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SUMMARY 


An integrated thermal-structural finite element approach for effi- 
cient coupling of thermal and structural analysis is presented. New 
thermal finite elements which yield exact nodal and element tempera- 
tures for one-dimensional ) inear steady-state heat transfer problems 
are developed. A nodeless variable formulation is used to established 
improved thermal finite elements for one-dimensional nonlinear transi- 
ent and two-dimensionai linear transient heat transfer problems. The 
thermal finite elements provide detailed temperature distributions 
without using additional element nodes and permit a common discretiza- 
tion with lower-order congruent structural finite elements. The accu- 
racy of the integrated approach is evaluated by comparisons with anal- 
ytical solutions and conventional finite element thermal- structural 
analyses for a number of academic and more realistic problems. Results 
indicate that the approach provides a significant improvement in the 
accuracy and efficiency of thermal-stress analysis for structures with 
complex temperature distributionsw 
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Chapter I 


^ INTRODUCTION 

The finite elecwmt method is one of the most significant develop- 
ments for solving problems of continuum mechanics. It was first 
applied by Turner el Ji* Cl)* in 1956 the analysis of complex 
aerospace structures. With increasing availability of digital 
computers) the method has become widespread and well recognized as 
applicable to a variety of continuum problems. Applications of the 
method to thermal problems were Introduced in the middle of 1960’s 
for the solution of steady-state conduction heat transfer [2l, 
Thereafter) extensions of the method were made to both transient and 
nonlinear analyses where nonlinearities may arise from temperature 
dependent material properties and nonlinear boundary conditions. 
Important publications of finite clement heat transfer analysis 
appear in references [3-12] , With these developments and consider- 
able effort contributed during the past decade ) the method has 
gradually increased in thermal analysis capability and become a 
practical technique for analyzing realistic thermal problems. 


*The numbers in brackets Indicate references. 
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1,1 Current Status of Thermal-Structural Analysis 

Thermal stresses induced by aerodynamic heating on advanced 

space transportation vehicles are an Important concern in structural 

design. Nonuniform heating may have a significant effect on the 
« 

performance of the structures and efficient techniques for determining 
thermal stresses are required. Frequently, the thermal analysis of 
the structure is performed by the finite difference method. 
Production-type finite difference programs such as MITAS and SINDA 
have demonstrated excellent capabilities for analyzing comrlex 
structures [13] . In structural analysis, however, the finite element 
method is favorable due to better capabilities in modeling complex 
structural geometries and handling various types of boundary condi- 
tions. To perform coupled thermal-structural analysis with efficiency, 
a computer program which includes both thermal and structural analysis 
codes is preferred, and a single numerical method is desirable to 
eliminate the tedious and perhaps expensive Cask of transferring 
data between different analytical models. 

Currently, the capabilities and efficiency of the finite element 
method is analyzing typical heat transfer problems such as combined 
conduction-forced convection is about the same as using the finite 
difference method [14]. With the wide acceptance of the finite 
element method in structures and its rapid growth in thermal analysis, 
it is particularly well-suited for coupled thermal-structural 
analysis. At present, several finite element programs which include 
both thermal and structural analysis capabilities exist; e.g. 

NASTRAN, ANSYS, ADINA and SPAR are widely used. These programs use 


3 


a common data base for transferring temperatures computed from a 
thermal analysis processor to a structural analysis processor for 
determining displacements and stresses. With the use of a common 
finite element discretization, a significant reduction of effort in 
preparing data is achieved and errors that ma/ occur by manually 
transferring data between analyses is eliminated. 

1.2 Needs for Improving Finite Element Methodology 

Although the finite element method offers high potential for 
coupled thermal-structural analysis, further improvements of the 
method are needed. Quite often, the finite element thermal model 
requires a finer discretization than the structural model to compute 
the temperature distribution accurately. Detailed temperature 
distributions are necessary for the structural analysis to predict 
thermal stress distributions including critical stress locations 
accurately. Improvement of thermal finite elements is, therefore, 
required so that a common discretization between the two analytical 
models can be maintalxied. 

Another need for improving the method includes a capability of 
the thermal analysis to produce thermal loads required for the 
structural analysis directly. At present, typical thermal-structural 
finite element programs only transfer nodal temperatures computed 
from the thermal analysis to the stiMctural analysis. These nodal 
temperatures are generally inadequate because additional information, 
such as element temperature distributions and temperature gradients, 
may be required to compute thermal stress distributions correctly, 


These needs are important in Improvement o£ finite element 
coupled thermal-structural analysis capability. The use of improved 
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thermal finite elements can reduce model size and computational 
costs especially for analysis of complex aerospace vehicle structures. 
Improved thermal elements will also have a direct effect in increasing 
the structural analysis accuracy through improving the accuracy of 
thermal loads. 

To meet these requirements for improved thermal- structural 
analysis and to demonstrate benefits that can be achieved, this 
dissertation will develop an approach called Integrated finite 
element thermal-structural analysis. First, basic concepts of the 
integrated finite element thermal-structural formulation are intro- 
duced in Chapter 2. Finite elements which provide exact solutions 
to one-dimensional linear steady-state thermal-structural problems 
are developed in Chapter 3. Chapter 4 demonstrates the use of these 
finite elements for linear transient analysis. Next, in Chapter 5 
a generalized approach for improved finite elements is established 
and its efficiency is demonstrated through thermal-structural 
analysis with radiation heat transfer. Finally, in Chapter 6 
extension of the approach to two dimensions is made with a new two- 
dimensional finite element. In each chapter, benefits of utilizing 
the improved finite elements are demonstrated by' both academic and 
realistic thermal-structural problems. 

Throughout the development of the improved finite elements, 
detailed analytical and finite element formulations are presented. 

Such details are provided in the form of equations, finite element 
matrices in tables and computer subroutines in appendices. 


Chapter 2 


AN INTEGRATED THERMAL-STRUCTURAL FINITE 
ELEMENT FORMULATION 

2.1 Basic Concepts 

Before applying the finite element method to thermal-structural 
analysis, It Is appropriate to establish basic concepts and procedures 
of the method. Briefly described, the finite element method is a 
numerical analysis technique for obtaining approximate solutions to 
problems by idealizing the continuum model as a finite number of 
discrete regions called elements. These elements are connected at 
points called nodes where normally the dependent variables such as 
temperature and displacements are determined. Numerical computations 
for each individual element generate element matrices which are then 
assembled to form a set of linear algebraic equations (for 
steady state problems) to represent the entire problem. These 
algebraic equations are solved simultaneously for the unknown 
dependent variables. Usually the more elements used, the greater 
the accuracy of the results. Accuracy, however, can be affected by 
factors such as the type of element selected to reprcssnt the con- 
tinuum, and the sophistication of element interpolation functions. 
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2.2 Element Interpolation Functions 

The first step after replacing the continuum model by a 
discrete number of finite elements Is to determine a functional 
relationship between the dependent variable within the element and 
the nodal variables. The function that represents the variation of 
a dependent variable Is called the Interpolation function. In thermal 
analysis, the element temperature T(x,.y,z,t) are generally expressed 
In the form 

T(x,y,z,t) - LNT^3t,y,z)J (T(t)} (2.1) 

where LN^(x,y,z)J denotes a row matrix of the element temperature 
Interpolation functions, and {T(t)} denotes a vector of nodal 
temperatures. Similarly, in a structural analysis, the element 
displacements, {6} , are expressed as, 

. {6(x,y,z,t)} » [Ng(x,y,z)] {5(t)> (2.2) 

where [Ng(x,y,z)] denotes a matrix of structural displacement inter- 
polation functions, and {6(t)} denotes a vector of nodal 
displacements. 

Usually, polynomials are selected as element interpolation 
functions and the degree of the polynomial chosen depends on the 
number of nodes assigned to the element. Regardless of the algebraic 
form, these interpolation functions have a value of unity at the node 
to which it pertains and a value of zero at other nodes. For example, 
linear temperature variation for a two-node one-dimensional rod 
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element vrith nodal temperatures and T 2 at x *0 (node 1) and 

X "L (node 2), respectively, can be written in the form 


T(x,t) - U - f 



Tj_(t) 


T^Ct) 


By comparing this equation with the general form of the element 
temperature variation, Eq. (2.1), the element interpolation functions 
are 


N^(x) ■ 1 - and N 2 (J^) " ^ 


These element interpolation functions, therefore, have the properties 
of N,( » 1 at node i and “ 0 at the other node. 


2.3 Finite Element Thermal Analysis 

Once the type of elements and their interpolation functions have 
been selected, the matrix equations expressing the properties of the 
individual element are evaluated. In thermal analysis, the method of 
weighted residuals [iS] is frequently employed starting from the 
governing differential equations. For condution heat transfer in a 
three-dimensional anisotropic solid 0 bounded by surface P 
(Fig. 1), an energy balance on a small element is given by, 


(-^ + 

dX 


3y 


3z 


■) + Q(x,y,z,t) » pc 


_^(x,y,z,t) 

3t 


(2.3) 


where q , q , q are components of the heat flow race per unit area, 
y z 

Q is the internal heat generation rate per unit volume, p is the 
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density, and c is Che specific heat. Using Fourier's Law, Che 
components of heat flow rate for an anisotropic medium can be written 
in the matrix form. 


*Ix 


^11 *' l 2 

*^13 


r *\ 

97/ax 

V 

% 

S J 

> m - 

*'21 * *^22 
*^31 ^32 

*^23 

^33 

•< 

9T/9y > 
3T/32 


where is the symmetric conductivity tensor. Figure 1 shows 

several types of boundary conditions frequently encountered in Che 
analysis. These boundary conditions are (1) specified surface 
temperatures, (2) surface heating, (3) surface convection, and 
(4) surface radiation: 


T - 


on S, 


qn +qn +qn «-q 
X X y y z 2 ^s 


on S, 


qn +qn +qn “ h(To - T„) on S- 
X X y y ^2 2 s 3 


qn + q.n +qn » ocT - aq„ on S, 
^x X T y 2 2 s 4 


(2.5a) 
(2.5b) 
(2.5c) 
. . (2.5d) 


where Tg is the specified surface temperature; n^, n^, n^ are the 
direction cosines of the outward normal to the sqrface, qg is the 
surface heating rate unit area, h is the convection coefficient, 

is the convective medium temperature, a is the Stefan-Boltsmann 
constant, e is the surface emissivity, a is the surface absorp- 
tivity, and .qj. is the incident radiant heat flow rate per unit 


area. 
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To apply Che finite element technique, the domain is first 
discretized into a number of elements. For an element with r 
nodes, the element temperature, Eq, (2,1), can be written in the 
form 

r 

T(x,y,z,t) - Z N. (x,y,z) T,(t) ’ (2.7a) 

i-1 

and the temperature gradients within each element are 


9T(x,y.z,t) . ' 

3x 3x 


(2.7b) 


3T(x,y,z,t) . J Ti(t) 


3y 


i-1 


3y 


(2.7c) 


3z i-1 3z 


(2.7d) 


These element temperature gradients can be wristen in the matrix 
form, 


3T (x,y,z,t) 
3x 



(x,y,z,t) 




(B(x,y,z)] (T(t)} (2.8) 


3T (x,y,z,t) 

3z 

where [B(x,y,z)] is the temperature-gradient interpolation 


matrix; 
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[B(x,y,z)] *• 


3Nj_ 

awo 

,iim 1 , _ ^ • • • 1 

3Nr 

1 • i • 

" ax’ 

ax 

ax 

3Ni 

awp 

T < 

3N^ 

1 . » . 

ay 

ay 

ay 

aN, 

aN, 

mmmmJSm • • * 1 

aN- 

az 

92 

as 


(2.9) 


and, therefore, the cowponenta of heat flow rate, Eq. (2.4), become 




i. V 




-[k3[Bj{T} 


( 2 . 10 ) 


where [k] denotes the thermal conductivity matrix. 

In the derivation of the element equations, the method of 
weighted residuals is applied to the energy equation, Eq. (2,3), for 
each individual element (e) , This method requires 


J 


n 


3x 3y ^ 3z 


■■ Q + pc 


11 

3t 


) N.dn 

i 


1,2 


( 2 . 11 ) 


After the integrations are performed on the first three terras by 
using Gauss’s Theorem, a surface integral of the heat flow across 
the eleraent boundary, T , is introduced, and the above equations 
becorae 


( [ (q -h) N. dr 



3N^ 9Nj_ 3N^ 
3x 3y 3z 



dn 


) 


q 


z 
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iff J- 

.. 



dn + 



dfi - 0 


( 2 . 12 ) 


where q is the vecter o£ conduction hcnt flux across the element 
boundary and il is a unit vector normal to the boundary. The 
boundary conditions as shown in Eqs. (2.5a -2.5d) are then imposed, 


^ I (t -n) dr - I qg + J \ dT 


+ J (creT^ - aq^) dF - J 


(e) 


9N^ 3N,j_ 9N^ 

' 9x' 9y 9z 


/• s 
^X 


> dn 


^ / 


r 

QN, dfi + pc N, dn 

„Ce) ‘ J (e) “ 

» n 


(2.13) 


By substituting the vector of heat flow rate, Eq. (2.10), the above 
element equations finally result in the matrix form, 


[C]{T} + [[K^] + [K^] + [K^]]{T} 

- {RJ + {R.} + {Rq}' + {R^} + {R^} 


(2.14) 


where [c] is the element capacitance matrix; [K^] , [K^J and 

[k ] are element conductance matrices corresponding to conduction, 

V * 

convection and radiation, respectively. These matrices are expressed 


as follows: 
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Ic) 

J,(.) '“t> LNtJ “« 

(2.15a) 

- 

IjjCe) ^®T^ dfl 

(2.15b) 

[K^] - 

J h {N.J,} [n^J dr 

S3 

(2. 15c) 

[Ky]{T} - 

1 oeT^ (Nj) dr 

®4 

(2.1Sd) 

The right-hand side of the discretized equation (2.14) contains 
heat load vectors due to specified nodal temperatures, internal heat 
generation, specified surface heating, surface convection and surface 
radiation. These vectors are defined by 

(y ■ 

- J <q-n) (N^) dr 

(2.16a) 

iRq) - 

i„(.) 

(2.16b) 

(R,> ■ 

L -•s'V'"' 
^2 

(2.16c) 

(y - 

J h T„ { N.J,} dr 

®3 

(2.16d) 

(Rj) - 

^4 

(2.16e) 

where q is the vector of conduction heat flux across 

, ^ 

boundary that 


is required to maintain the specified nodal temperatures. 


2.4 Finite Element Structural Analysis 


In a finite element structural analysis, element matrices may 
bo derived by the method of weighted residuals, or by a variational 
method such as the principle of minimum potontiol energy tl7-19] . 

For simplicity in establishing these element matrices and understand- 
ing general derivations, the last approach is presented herein. 

The basic idea of this approach is to derive the static equilibrium 
equations and then include dynamic effects through the use of 
D'Alembert's principle. 

Consider an elastic body in a three-dimensional state of stress. 
The iMtemal strain energy of an element (e) can be written in a 

t 

form, 

U j'^(e)[e-eoJ io) dn (2.17) 

* 

where is the element volume, {cr} denotes a vector of stress 

components; [ej and [eqJ denote row matrices of total strain and 

initial strain components, respectively. Using the stress-strain 
relations, 


{a} - [d] (e - Eq} (2,18) 

where [d] is the elasticity matrix, the internal strain energy 


U - 


1 

2 



[e - EqJ [d] {e - Eq} dn 


becomes 



or 


IS 


“ ■ 1 1(e) W t®l M <ia - f UJ (Dl Ce,) dfl 

+ 1 Icpj [Dl (ej) <18 (2.19) 

For each elemont, the potential energy of the external forces 
may result from body forces and boundary surface tractions* The 
potential energy due to body forces can be written in a form, 

‘ Vq - - I UJ {f} dft (2*20) 

where (f) denotes a vector of body force components. Similarly, 
the potential energy due to surface tractions is, 

Vg - - UJ U) dr . (2.21) 

where {g} denotes a vector of surface traction components, and 
denotes the element boundary. The total element potential 
energy, tt^ , the sum of the internal strain energy and the potential 
energy of the e'xternal forces is, 

(e) ^®J “ L(e) 1-®^J 

+ 1 I„(e) '■>) '=0> “ - 4(s) «« 

- I (e) LOJ (8) 


( 2 . 22 ) 


OF POOH QUhlUy 



where u, v, w are components of displacement In the three coordinate 
directions. The vector of strain components can be computed from 

8u 
3x 

Iz 
ay 

aw 

3z 

3u ^ 3v 
9y 9x 

9v ■ 9w 
9z 9y 

9u ■ 9w 
9z 9x 

where [Bg] is the strain-displacement interpolation matrix. By 
substituting the element displacement vector, Eq. (2.23), and the 
vector of strain components, Eq. (2.24) into Eq. (2.22), the total 
element potential energy is expressed in tenns of the nodal displace- 
ment vector {6} as 
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OF poon QUAury 

% - I l«J f (0) '“si'' '“1 'M <'o> 

+ T j (e) L'oJ ■*“ - I (e) '*) 

n lill! 

- {6} |^(e) [Ns]"^ {8} dr (2.25) 

The principle of minimum potential energy requires, 

aiTe 

— ^ • 0 
3{6} 

which yields the element equilibrium equations, 

[Kg] {6} - {F^} + {Fg} + {Fg} + {Fj} (2.26) 

where [Kg] is the element stiffness matrix defined by 

[Kg] - [Bg]^ [D] [Bg] dfi (2.27a) 

The right hand side of the equilibrium equations contains force 
vectors due to concentrated forces, body forces, surface tractions 
and initial strain, respectively. The nodal force vectors due to 
body forces and surface tractions are 

(2.27b) 
(2.27c) 


^^B> “ Le) tf> dn 

fFs> “ J (e) {8} dr 
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For initial strains from thermal effects, the corresponding nodal 
vector (F^} is due to the change of temperature from a reference 
temperature of the zero-stress state and may bo written as 

" J (e) (2.27d) 

n I 

where {a} is a vector of thermal expansion coefficients, T is 
the element temperature distribution, and is the reference 

temperature for zero stress. 

For elastic bodies subjected to dynamic loads, the effects of 
inertia and damping forces roust be taken into account. Usir.g 
D’Alembert's principle, the inertia force can be treated as a body 
force given by 

{f} - - p{6} (2.28a) 

where p is the mass per unit volume. By using element displacement 
variations, Eq. (2.23), this inertia force is expressed in terms of 
nodal displacements *^3 

{f} - - P [N J it} (2.28b) 

S t 

Similarly, the damping force which is usually assumed to be propor- 
tional to the velocity can be expressed in the form, 

{f} - - p CNg] (2.28c) 

where p is a damping coefficient. By substituting these inertia 
and damping forces, Eqs. (2.28b -2.28c), into Eq. (2.27b), the equi- 
valent nodal body forces shown in Eq. (2.27b) become 
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■ I (4) ““ - 1(e) 1N,1 d!! (5) 

Q H 

- j (e) [Na]'’^p [Ng] dfl {6} (2.29) 

Finally, by using the static equilibrium equations, Eq. (2.26), with 
the above equivalent nodal body force, the basic equations of struc- 
tural dynamics can be written in the form, 

[M] {6} + (Cg] {S} + [Kg] is} - (Fc) + {Fg} + {Fg} + {F^} (2.30) 

where [m] and tCg] are the element, mass and damping matrices, 
respectively, and defined by 


[M] - ’ 

J 

[(e) P 

[Ng] dJ2 

(2.31a) 

[Cs] - 

w 

L(e) " 

[Ng] dn 

(2.31b) 


In a general formulation of transient thermal-stress problem, 
the heat conduction equation (2.3) contains a mechanical coupling 
term in addition [16], This coupling term represents the mechanical 
energy associated with deformation of the continuum and in some 
highly specialized problems (see Ref. 16) can affect the temperature 
solution. In most of engineering applications, fortunately, this term 
is insignificant and is usually disregarded in the heat conduction 
equation. This simplification permits transient thermal solutions 
and dynamic structural responses to be computed independently. 

For a structural analysis where the inertia and damping effects 
are negligible, the static structural response, Eq. (2.26), can be 
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computed at selected times corresponding to the transient thermal 
solutions. Such a sequence of computations, widely used in thermal- 
structural applications, is called a quasi-static analysis. Results 
of temperatures directly enter the structural analysis through the 
computation of the thermal nodal force vector, Eq. (2.27d). Tempera- 

I 

tures also have an indirect effect on the analysis through the 
structural material properties, since the elasticity matrix [D] 
and the thermal expansion coefficient vector {a} are, in general, 
temperature dependent. Temperature dependent properties may result 
in a variation of the structural element stiffness matrix, Eq. (2.27a), 
throughout the transient response. 

2.5 Integrated Approach 

The representation of the element temperature distribution in 
the computation of structural nodal forces is an Important step in 
the coupled thermal-structural finite element analysis. In typical 
production-type finite element programs, element nodal temperatures 
are the only informatioii transferred from the thermal analysis to 
the structural analysis. This general procedure is shown schemati- 
cally in Fig. 2(a) and herein is called the conventional finite 
element approach. Since the conventional thermal analysis only 
provides nodal temperatures, an approximate temperature distribution 
is assumed in the structural analysis which results in a reduction 
in accuracy of displacements and thermal stresses. 

To improve the capabilities and efficiency of the finite 
element method, an approach called integrated thermal-structural 
analysis is developed as illustrated by Fig. 2(b). The goals of 


CONVENTIONAL INTEGRATED 
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Fig. 2. Conventional versus Integrated thermal and 
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the integrated approach are to; (1) provide thermal elements which 
predict detailed temperature variations accurately, (2) maintain 
the same discretization for both thermal and structural models with 
fully compatible thermal and structural elements, and (3) provide 
accurate thermal loads to the structural analysis to improve the 
accuracy of displacements and stresses. 

These goals of the Integrated approach require developing new 
thermal finite elements that can provide higher accuracy and effi- 
ciency than conventional finite elements. The basic restriction on 
these new thermal elements is the required compatabllity with the 
structural elements to preserve a common discretization. Detailed 
temperature distributions resulting from the improved thermal finite 
elements can provide accurate thermal loads required for the 
structural analysis by rigorously evaluating the thermal load 
integral, Eq. (2.27d). 
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Chapter 3 

EXACT FINITE ELEMENTS FOR ONE-DIMENSIONAL 
LINEAR THER^^AL-STRUCTURAL ..EMS 

In general, polynomials are selected as element Interpolation 
functions to describe variations of the dependent variable within 
elements. In one-dlmenslonal analysis, the simplest polynomial which 
provides a linear variation within an element Is of the first order, 

- Cl + C 2 X (3.1) 

where <j> denotes the dependent variable such as temperature or 
displacement; Ci and C 2 denote constants, and x is the coor- 
dinate of a point within the element. A finite element with two 
nodes is formulated by imposing the conditions at nodes, 

(j)(x =0) = <t>i (fi(x=L) » ^2 (3.2) 

where L is the element length; <|)i and <(>2 are nodal values at 
node 1 and 2, respectively. The dependent variable, therefore, can 
be written in terms of nodal values as 

<}> = (1 - f ) <(>1 + (^) < l >2 

or in the matrix form, 
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- [nJ C<j>} (3.3) 

where [nJ is the row matrix of element interpolation functions. 

The type of finite element where the dependent variable is assumed 
to vary linearly between the two element nodes is often used in one- 
dimensional problems and is called a conventional finite element 
herein. Vflth the linear approximation, a large number of elements 
are required to represent a sharply varying dependent variable. In 
some special cases, however, conventional finite elements can provide 
exact solutions when the solutions to problems are in the form of a 
linear variation. For example, a linear temperature variation is the 
exact solution of one-dimensional steady-state heat conduction in a 
slab; therefore, the use of the conventional finite element leads to 
an exact solution. Further observation [20] has shown that, under 
some conditions, exact nodal values are obtained through the use of 
this element type. Temperatures for steady-state heat conduction 
with internal heat generation in a slab and deformations of a bar 
loaded by its own weight are examples of this case. In the past, 
the capability of conventional finite elements to provide exact 
solutions has been regarded as a property of the particular equation 
being solved and not applicable to general problems. 
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In thin chapter, finite elements that provide exact solutions 
to one-dimensional linear steady-state thermal-structural problems 
are given. The fundamental approach in developing exact finite 
elements is based on the use of exact solutions to one-dimensional 
problems governed by linear ordinary differential equations. A 
general formulation of the exact finite element is first derived and 
applications are made to various thermal-structural problems. 
Benefits of utilizing the exact finite elements are demonstrated 
by comparison \rLth results from conventional finite elements and 
exact solutions. 


3.1 Exact Element Formulation 

In this section, a general derivation of exact finite elements 
is given. Exact finite elements for various thermal and structural 
problems are derived and described in detail in the subsequent 
sections. Consider an ordinary, linear, nonhomo geneous differential 
equation, 


a illi + a 

“ dx" dx^“l 


+ + “o* ■ 


(3.4) 


where x is the independent variable, (|)(x) is the dependent 
variable, a^, i»0, n are constant coefficients, and r(x) is 

the forcing function. A general solution to the above differential 
equation has the form 

n 

<|>(x) - E C, f . (x) + g(x) 
i-1 


(3.5) 
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whore Cj aro arbitrary constanto, typical functions in 

the homogeneous solution and g(x) is a particular solution. For 
example, a typical one-dimensional steady-state thermal analysis is 
governed by second order differential equation of the form of 
£q* (3.4) and has a general solution 

(}i(x) •• CjL f^Cx) + C 2 f 2 (x) + g(x) (3.6) 

By comparing this general solution with the solution in the form of 
polynomials used to describe a linear variation of dependent variable 
in the conventional finite element, Eq. (3.1), basic differences 
between these two solutions are noted: (1) the function f;j^(x) in 

the general solution to a given differential equation can be forms 
Ocher than the polynomials, and (2) the general solution contains a 
particular solution g(x) which is known in general and depends on 
forcing function r(x) on the right hand side of the differential 
equation (3.4). 

3.1.1 Exact Element Interpolation Functions and Nodeless 
Parameters 

Once a general solution to a given differential is obtained, 
exact element interpolation functions can be derived. For a typical 
finite element with n degrees of freedom, n boundary conditions 
are required. With the general solution shown in equation (3.5), 
the required boundary conditions are 

<j)(x^) » i “ 1, 2, , n (3.7) 

where x. is the nodal coordinate and $ . is the element nodal 
i X 


unknown at node 1. After applying the boundary conditions, the 
exact element variation of ‘{'(x) has the form, 
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n 

(x) - G(x) + E N.h 
i»l ^ ^ 


where H^(x) I's the clamant interpolation function carrosponding 
to node 1, The function G(x) is a known function associated with 
the particular solution. In general, this function can be expressed 
as a product of a spatial function Nq(x) and a scalar term 
which containa a physical forcing parameter such as body force, 

I 

surface heating, etc.; 

G(x) - NqCx) <^0 

and, therefore, the exact element (j>(x) variation becomes. 


i|)(x) 




(3.8a) 


or in the matrix form 


4>(x) - [Nq ^1^2 


LNq 


'J’O 

T"^ 


(3.8b) 


T" J 

Note that the element interpolation function has a value 

of unity at node i to satisfy the boundary conditions, Eq. (3.7), 
thus the spatial function Nq(x) must vanish at nodes, Since the 
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com <|)q io a knovm quanclcy and neither tolatos to the olcmont 
nodal coordinatoo nor io idontifiod with the element nodeo, it is 
called a nodeloso parameter! Llkowise, the corresponding spatial 

o 

function Nq(x) Is called a nodelcss interpolation function. 

Comparison between element variations of a typical nodelcss para- 

* 

meter finite element, Eq. U.S), and the conventional linear finite 
element) Eq. (3.3), is shown in Fig. 3, 


3.1.2 Exact Element Matrices 

After exact element interpolation functions are obtained, the 
corresponding element matrices can be formulated. For the governing 
ordinary differential equation, Eq. (3.4), typical element matrices 
can be derived (see section 2.2) and element eqvmtions can be written 
in the form. 
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> * 


(3.9) 


where K^, j , i, j »0, n are typical terms in the element stiffness 
matrix; F^, 1 "0, n are typical terms in the element load vector, 
ij)^, i “1, n are the element nodal unknowns, and (^q is the element 
nodeless parameter. Since the element nodeless parameter is known, 
the above element equations reduce to 
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X 


Fig. 3.^ Compgalson of conventional and nodeless parameter 
elements , 
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^nO 
> • 


(3.10) 


• 3.2 Exact Finite Elements in Thermal Problems 

In one-dimensional linear steady-state thermal problems, typical 
governing differential equations can be derived from a heat balance 
on a small segment in the form, 

S f) + *!<=<> S + “oW 

where T denotes the temperature, x denotes a typical one- 
dimensional space coordinate in cartesian, cylindrical or spherical 
coordinates; a^, i “0, 1,2 are variable coefficients, and r(x) 
is a function associated with a heat load for a given problem. A 
general solution to the above differential equation has the form, 

T(x) = f^(x) + + g(x) (3.12) 

where f^(x) and linearly independent solutions of the 

homogeneous equation, and C 2 are constants of integration, 

and g(x) is a particular solution. Since the particular solution 
g(x) is known, the above general solution has two unknowns to be 
determined. A finite element with two nodes, therefore, can be 


formulated using the conditions. 
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T(Xj^) - (3,13a) 

T(x2) • T 2 (3.13b) 

where x^, i»l> 2 are nodal coordinates and i "• 1, 2 are the 

nodal temperatures. Imposing these conditiohs on the. general solu- 
tion yields two equations for evaluating and C 2 , 

T(Xi) " + ^2 ^2^*1^ ®^^1^ 

TCX 2 ) " T 2 - ^2 ^2^’*^2^ ®^^2^ 

or in matrix form 


f 'I t- ^ 



to 



i 

Ti - g(Xj^) 

‘ r 




^2 


T 2 -- g(x2) 




After and C 2 are determined and substituted into the general 

solution, Eq. (3.12), the exact element temperature variation can 
be written as 


T(x) - Nq(x) Tq + Nj^(x) + N 2 (x) T 2 (3.14a) 


or in the matrix form. 


T(x) - LNq N 2 J ^ 




"2 


(3.14b) 


where Nq(x) is the nodeless interpolation function and Tq is the 
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tiodeless parameter; Nj^(x) and N 2 (x) are element interpolation 
functions corresponding to node 1 and 2, respectively, these element 
interpolation functions including the nodeless parameter are known 
functions defined by 


f. (x) f,(x 2 ) - 

N, (x) “ — ■= =--^ 

^ W 


(3,15a) 


f-(x-) f 2 (x) - fj_(x) 

N, (x) — 

^ W 


(3.15b) 


Nq(x) 


f,(x.) g(x„) - f,(x 2 ) g(x. ) 
g(x) ^ ^ ^ f,(x) 

w 


^1^^2^ g(Xi) - f^(x^) g(x2) ^ ^ ^ 
W 


(3.15c) 


where W = ^2^^2^ ” ^1^’^2^ ’ 


Using the exact element interpolation functions shown in 
Eq. (3.14), and the governing differential equation, Eq. (3.11), 
element matrices can be derived through the use of the method of 
weighted residuals; 


P 2 


X, 


Performing an integration by parts on the first term and substituting 
for element temperature in teiniis of the interpolation functions, 

Eq. (3.14), yields element equations in the form. 
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[[K^] + [Ky] + [Khl] 


r ^ 





V. J 


{Qg} + {Q} 


(3.17) 


where » and [K^] are trhe element conductance matrices 

associated with the second, the first, and the zero-order derivative 
term on the left hand side of the governing differential equation 
(3.11), respectively; {Q^} is the element vector of conduction 
heat flux across element boundary, and {Q} is the element load 
vector from the heat load r(x’) in the governing differential 
equation. These matrices are defined as follows: 


!k,1 


f dNi I dN| 
^2 'te* 


(3.18a) 


[K ] 

V 


^1 


(3.18b) 


i\] “ \ {N} [nJ dx 


(3.18c) 


X.. 


'V = it- " 1 


(3.18d) 
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{Q} 




(3.18e) 


Depending on the complexity of element interpolation functions, the 
element matrices may be evaluated in closed fotm or they may require 
numerical integration. However, after the element matrices are 
computed, typical element equations can be written in the form, 



(3.19) 


Since the nodeless parameter is known, the first equation is uncoupled 
from the nodal unknowns in the second and third equations. Thus, 
the exact element matrices have the same size as of the conventional 
linear finite element and element equations can be written as 



Note that, in general, the above conductance matrix is an 
asymmetric matrix. This asymmetry is caused by the conductance 
matrix [k^] shown in Eq. (3.18b) associated with the first-order 
derivative in the governing differential equation, Eq. (3.11). To 
obtain a symmetrical conductance matrix, the first-order derivative 
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is eliminated by casting the governing differential equation in self- 
adjoint form, 


^ [P(x) + Q(x) T - R(x) (3.21) 

where 

P(x) - exp(J[(aj^ ■**. ”dx^ ^ (3.22a) 

(3.22b) 

(3.22c) 


Q(x) 


®2 


R(x) - ^ 


Element matrices can then be derived using the method of weighted 
residuals in the same manner as previously described. In this case, 
element equations have the form. 


[[KeJ + [%]] 


> » {Qc> + tQ> 


(3,23) 


where the conductance matrices and heat load vectors are defined 
by 



(3.24a) 
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n^ P'lOR OUAUTY 


[Kj^] - Q {N} [nJ dx 


'«c> 


< 

dx 


X, 


X, 


{Q} - 1 R {N} dx 


(3.24b) 


(3.24c) 


(3.24d) 


Similarly, element equations for the two nodal unknowns are 



^11 

7^1 

K) 

I 


^12 

^22 

where , i,j 

= 0,1, 

coefficients in 

the c 


f - 

*^1 
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> - T ^ 

\o' 

> 

^2 

u 

^20 
^ - 


(3.25) 


Q Nj dx (3.26) 

i,j=0,l,2 

An additional advantage of using the self-adjoint differential 
equation is that the coefficients and K 2 q shown on the right 

hand side of Eq. (3.25) are identically zero. This result can be 
proved by observing that the element interpolation function N^, 


x„ 


dN. dN. f 

dx + 

dx dx 


K, 


ij 


- P 
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i »1,2 are the solution of the homogeneous differential equation, 
Eq. (3.21), because is a linear combination of the homogeneous 

solutions f^(x) and f 2 (x) as shown in Eqs. (3.15a-b), i.e, 

^[P^]+Q»l-0 1-1.2 


Multiplying this equation by the nodeless parameter interpolation 
function Nq and performing integration by parts on the first term 
yields 





dN. dN 

p I ■■■ I II «—■ — ft H V 

dx dx ^ ^ 


Q Nq dx 


X, 


Then since the nodeless interpolation function Nq vanishes at 
nodes, i.e. at the coordinates x^ and X 2 » the above equations 
yield 


K^q - 0 i -1,2 

and the element equations, Eq. (3.25), become 


1 

1 

1 

|-» 

to 

I 


f 

"1 

> » {Q^} + ^ 

r - 

^1 

> 

^12 

K22 


."2 


^2 

•• ^ 


After element nodal temperatures are computed, exact temperatures 
within an element can be obtained using the exact element temperature 
variation, Eq. (3.14). 
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To demonstrate the exact finite element formulation previously 

derived, exact finite elements for eight heat transfer oases in 

several solids of different shapes and a flow passage (Fig, 4) are 

presented. In the first seven oases, heat transfer may consist of; 

(1) pure conduction, (2) conduction with internal heat generation, 

» 

(3) conduction with surface heating, and (4) conduction with surface 
convection. Case eight is a one-dimensional flow where heat transfer 
may consist of fluid conduction and mass transport convection with 
surface heating or surface convection. For these cases, the boundary 
conditions considered are: 



T “ Constant 

(3.23a) 

or 

- q 

dx ^ 

(3.28b) 

or 

- k g - h(T -TJ 

(3.28c) 


where k is the material thermal conductivity, q is the specified 
surface heating rate per unit area, h is the convection coefficient, 
and T„ is the convection medium temperature. In each case, the 
derivation of exact finite elements for appropriate heat transfer 
cases are given for clarity. Governing differential equations and 
the corresponding nodeless parameters, exact element interpolation 
functions, and element matrices for all cases are shown in Tables 1 
and 2 and Appendices A and E. 

3.2.1 Rod and Slab 

A rod element with arbitrary cross-sectional area A, circum- 
fer... -ial perimeter p and length L as shown in (Fig. 4, Case 1) 





Governing Self-Adjoint Differential Equations 
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*Corabined conduction and mass transport convection where P = exp(-mcx/kA) 
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Table 2 



Nodeless Parameters 
* 

for Thermal Problems 




^0 


Case 

Convection 

Source 

Surface Flux 


(b) 

(o) 

(d) 

1 

T 

2b! 

qpL^ 

JU 

^00 

2k 

2kA 

2 


fib! 




2k 


3 


fi^ 




4kw 


4 






6k 


5 

T 

fib?. 

t2 

3L. 


■**09 

2k 

2kt 

A 

T 

fib! 

u2 

-aL- 

u 

■^oa 

4kw 

4ktw 

7 


fia! 

fia! 



k 

kt 

8 

T 


S£b 


^00 


ibc 


where 


w » ln(b/a) 
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is subjected to internal heat generation, surface heating, and 
surface convection. Governing differential equations for each heat 
transfer case in self-adjoint form are shown in Table 1, For example, 
the governing differential equation for the case of conduction with 
surface convection is 


- (kAg) + hpT - hpT„ (3.29) 

where Ic is the material thermal conductivity, h is the convection 
coefficient, and is the convective medium temperature. A general 

solution to the above differential equation is 

T(x) ■ sinh mx + C 2 cosh mx + 

where m ■ /hp/kA, and and C 2 are unknown constants. Applying 

the boundary conditions at the nodes, 

ff 

T(x“0) - and T(x-L) ■ T 2 

the two unknown constants are evaluated and the above solution 
becomes 


T(x) 


(1 - 


sinh m(L-x) 
sinh mL 


sinh mx .. _ 
sinh mL'^ 


00 


. ,sinh m(L-x) > „ , ,slnh mx > 

sinh mL ' '^1 ^sinh mL'' 2 


(3.30) 


This exact element temperature variation can be written in the form 
of Eq. (3.14) where the element interpoation L ;nctlons and the node- 
less parameter are; 
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No<x) 


/1 Blnh m(L-x) slnh mx ^. _ ^ * 

“ slnh raL " sinh mL*^’ 0 




alnh m(L~x) , 
slnh mL ’ 


N2(x) 


slnh nuc 
slnh mL 


(3.31) 


As described In the previous section, element equations for a 
typical self-adjoint differential equation have the form of Eq. (3.23), 
and usln({ the definitions of the element matrices shown In Eq. (3.24), 
the element matrices for this problem are: 


L 

ty . f kA<f> L^J d* (3.32a) 

■'o 

[Kjj] - 1 hp {N} [Nj dx 

•^0 ’ 

L 

(Q) - f hpT„ {N} dx 
'*0 

where [k 1 and [iCu] are conductance matrices corresponding to 

C » 

conduction and convection, respectively, and {Q} Is the load vector 
due to surface convection. With the exact Interpolation functions 
shown in Eq. (3.31), the above element matrices can be evaluated in 
closed form. Exact nodal temperatures and element temperature varia- 
tion can then be computed using Eqs. (3.27) and (3.30), respectively. 

For the cases where the rod is subjected to an internal heat 
generation or specified surface heating, exact element interpolation 
functions and element matrices can be derived in the same manner as 


(3.32b) 


(3.32c) 
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described above. It should be noted that only the conductance matrix 
associated with conduction and heat load vectors corresponding to 
Internal heat generation or surface heating exist In the two cases. 
The exact conductance matrix and heat load vectors are found to be 
Identical to those from the conventional linear element, Therefore, 
exact nodal temperatures can also be obtained through the use of 
the conventional linear finite element In such cases. However, since 
the linear temperature variation Is not an exact solution to these 

problems, the conventional linear finite element can not provide the 

« 

exact temperature distribution within the element. 

The derivation of exact finite elements for one-dimensional heat 
transfer in a slab follows the derivation for the exact rod element. 

A slab with thickness L subjected to an internal heat generation 
(Fig. 4, Case 2) where both sides of slab may be subjected to a 
specified surface heating or surface convection. In Table 1, the 
governing differential equations are shown only for the case of pure 
conduction and conduction with internal heat generation because the 
effects of surface heating and surface convection enter the problem 
through the boundary conditions. For example, a governing differen- 
tial equation describing heat conduction in a slab with specified 
temperature T^ at x = 0 (node 1) and surface convection at x = L 
(node 2) is 




0 


(3.33) 


where k denotes the material thermal conductivity. After solving 
for the general solution to the governing differential equation above 



and applying nodal temperatures as boundary conditions at x ■ 0 
and X - L, the exact element temperature variation is (see 
Appendix A) 
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T(x) - (1 - Tj + (|) T 2 - LNj_(x) N 2 (x)J 



(3.34) 


With the corresponding element conductance matrix shown in Appendix B, 
exact element equations for this problem are 



since at x ■ L (node 2) the boundary condition is 


^ (x-L) 

dx 


h(T2 -T„) 


where h is the convection c ©'Efficient and T„ is the surrounding 
medium temperature, therefore, the above element equations become 



the exact nodal unknown T 2 can then be computed and the exact 
element tempterature distribution is obtained using equation (3.34). 

The same procedure can be applied for the case when the slab is 
subjected to surface heating. In this case, the boundary condition 


is 
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where q denotes the specified surface heating, \\fhen the slab 
consists several layers with different thermal conductivities, an 
exact element can be used to represent each layer. If the slab is 
subjected to surface heating or surface convection in addition, the 
above procedure applies for the elements located at the outer 
surfaces. 

3.2.2 Hollow Cylinder and Sphere 

A thermal model of a hollow cylinder with radial heat conduction 
subjected to an internal heat generation is shown in Fig. 4, Case 3. 
Specified heating or surface convection are considered through the 
boundary conditions at the inner and outer surfaces of radii a and 
b, respectively. Governing differential equations corresponding to 
each heat transfer case are provided in Table 1. For example, the 
governing differential equation for the case of pure conduction is 

k^[rgj..O (3.36) 

where k is the material thermal conductivity, and r is the radial 
coordinate. A general solution to the above differential equation is 

T(r) = + C 2 In r 

Nodal temperatures are imposed on the element boundary conditions, 

T(r =»a) = and T(r =b) = T 2 
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and the exact element variation is obtained as (see Appendix A), 

Tfr) « I ln(b/r) ln(r/a) i 

where w ■ ln(b/a). Note that the exact element variation for this 
case is completely different from the linear element variation, there' 
fore, the conventional linear finite element can not provide exact 
element or nodal temperatures. Applying the method of weighted 
residuals to the governing differential equation, element equations 
are 


j 


(3.37) 


[Kj.] {T} » {Qj,} 


(3.38a) 


where 



(3.38b) 


(3.38c) 


Using the exact element interpolation functions shown in Eq. (3.37), 
element equations for this case are 



dr 


kb 


dr 


(3.39) 
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When the cylinder is subjected to surface heating or surface 

> 

convection, the same procedure previously described for the slab 
can be used. For example, in case of convection heat transfer on 
the outer surface, the boundary condition is 

r - b; -k||-h(T 2 -Tj 

where h is the convection coefficient and T„ is the surrounding 
medium temperature. Thus, the element equations, Eq. (3.39), become 

(3.40) 

« 

Exact finite elements can be formulated for conduction heat 
transfer in the radial direction of a hollow sphere with internal 
heat generation. A thermal model of a hollow sphere with inner and 
outer surface radii a and b, respectively, is illustrated in 
Fig. 4, Case 4. The hollow sphere may be subjected to surface 
heating or surface convection on both inner and outer surfaces. 

For heat conduction with internal heat generation, the governing 
differential equation is (see Table 1) 

(3.«) 

where k is the material thermal conductivity, Q is the heat 
generation rate per unit volume, and r is the independent variable 
representing the radial coordinate. A general solution to this 



differential equation is 
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Or^ 

T(r) - "7 + “ V 

Due to the presence of the particular solution in the above general 
solution, a nodeless parameter exists, and the exact element 
variation is written in the form 

T(r) - NQ(r) Tq + Nj^(r) + N 2 (r) (3.42a) 

where the element interpolation functions including the nodeless 
parameter are: 


NgCr) - -^(r-a)(b-r)(r+a+b); 


M (j.) , aO^rl 
r(b-a) 




(3.42b) 


Element matrices can be derived using the method of weighted residuals 
and element equations are resulted in the form 


[K_] {T} = {Q.} + {Q} 


(3.43a) 


where these element matrices are defined by: 


t-Kj 


f I* dr 

dr ‘•dr’' 


(3.43b) 


{Qe> = ^ 


kr2f N 
dr 


(3.43c) 
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b 

{Q} - f Q{N} dr (3.43d) 

a 

If surface heating and surface convection are applied on the inner 
and outer surface, the same procedure described for the cylinder is 
required. 

3.2.3 Thin Shells 

Three thermal models of thin shells of revolution with cylin- 
drical, conical and spherical shapes are presented (see Fig. 4). 

These shells may be subjected to thermal loads such as surface 
heating, surface convection, and internal heat generation as shown 
in Fig. 4, Cases 5-7. In Case 5, a cylindrical shell of radius a, 
thickness t and meridional coordinate s is considered. Governing 
differential equations corresponding to different thermal loads are 
shown in Table 1. These governing differential equations are in the 
same form as for the rod element (Case 1). Therefore, the exact 
rod element interpolation functions and element matrices previously 
derived can be modified and used for the exact cylindrical shell 
element . 

A truncated conical shell element with thickness t is shown 
in Fig. 4, Case 6. Governing differential equations corresponding to 
internal heat generation and surface heating are given in Table 1. 
These differential equations are in the same form as for the hollow 
cylinder (Case 3) with surface heating, and therefore, element 
interpolation functions and element matrices are similars. For 
the case of the shell subjected to surface convection, a form of 
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nonhomogeneous modified Bessel's differential equation results, 


d^T I 1 dT h , _ h 

s ds kt kt " 


(3.44) 


A general solution to the above differential equation Includes 
modified Bessel functions of the first and second kind of order zero. 
A nodelesB parameter also exists in this case due to the nonhomo- 
geneous differential equation. Applying nodal temperatures as the 
boundary conditions at s ■ a and s » b, exact element interpo.la- 
tlon functions are obtained as shown in Appendix A. 

Fig. 4, Case 7 shows a truncated spherical shell with radius 
a and thickness t. The spherical shell may be subjected to 
internal heat generation or surface heating. Governing differential 
equations corresponding to these thermal loads are in the form of 
Legendre's differential equation of order zero. For example, the 
governing differential equation for the case of uniform surface 
heating q is 


(1 - 


2, d^T 

1 ) 

dn 


- 2n 


dn 


2 

_ S§_ 
kt 


(3.45) 


where n = sin (s/a) . A general solution to the above differential 
equation is 


T 





(3.46) 


where and are unknown constants. By imposing nodal 

temperatures as element boundary conditions at s = 0 and s = L, 
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exact elemGnt interpolation functions are obtained as shown in 
Appendix A. 

Due to the complexity of the exact element interpolation func- 
tions that arise from the truncated conical shell with surface 
convection and the truncated spherical shell, the corresponding 
element matrices in closed form are not provided. The element 
matrices, if desired, can be obtained using the element matrix 
formulation shown in Equations (3.14a-d) and performing the integra- 
tions numerically. 


3.2.4 Flow Passage 

A thermal model of fluid flow in a passage with conduction and 
mass transport convection is illustrated in Fig. 4, Case 8. The 
fluid may be heated by surface heating, or surface convection. 
Governing differential equations corresponding to these heat transfer 
cases are given in Table 1. For simplicity, consider the case without 
heat loads where the governing homogeneous differential equation is 
given by 


d r , , dT 1 , • 

- -T— [kA -T-] + me 
dx dx"* 


dx 


(3.47) 


where k is the fluid thermal conductivity, A is the flow cross- 
sectional area, m is the fluid mass flow rate, and c is the fluid 
specific heat. A general solution to this differential equation is 


T(x) “ + C 2 exp (2ax) 


where and C 2 are arbitrary constants and a =» mc/2l<A. An 

exact finite element with length L and nodal temperatures T^^ and 


S3 


at X - 0 and x ■ L, respectively, can be formulated, The 
exact element temperature variation is 


T(x) » U - “ 
1 




(3.48.a) 


As previously described, the appearance of the first-order derivative 
term in the governing differential equation results in an unsynnnetrical 
conductance matrix (see Eq, (3.18b)). In this case, the corresponding 
element conductance matrices are 


[K 1 - kAa 
c 


(e 


2aL 


(e 


2aL 


+ 1 ) 
- 1 ) 


1 



-1 




(3.48b) 


(3.48c) 


where tK^]‘ and [K^] denote conductance matrices representing 
fluid conduction and mass transport fluid convection, respectively. 

It has been shown that if the conventional finite element with 
an optinum upwind weighting function is used, exact temperatures at 
nodes can also be obtained [21 ]. With upwind weighting functions 
the element temperature variation is expressed as, 


T(x) 


[1 - f + F(x) 



(3.49a) 


where F(x) is the optimum upw"J.nd weighting function defined by 
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F(x) * [coth (2aL) 




With these element Interpolation functions, element conductance 
matrices corresponding to the fluid conduction and mass transport 
convection are 


upwind 


^ upwind 



+ ^ (coth (CLl) - M 
2 ^ 




(3.49b) 


(3.49c) 


It can be shown that the combination of these element conductance 
matrices are Identical to those obtained from the exact finite 
element, Eqs. C3.48b-c). Therefore, the conventional finite element 
with the optimum upwind weighting function provide exact nodal 
temperatures. However, since the upwind element temperature varia- 
tion differs from the exact element temperature variation shown in 
Eq. (3.48a), the finite element with the optimum upwind weighting 
function does not provide the exact temperature variation within an 
element . 


3.3 Exact Finite Elements in Thermal-Structural Problems 

With the general exact finite element formulation described 
in section 3.1, exact structural fxnite elements can be developed 
for problems governed by ordinary differential equations. For 
example, exact finite elements for a rod loaded by its own weight 
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or a beam with a dlstribuCcd load can be formulated. However, for 
the purpose of demons Crating benefits on exact finite elements in 
coupled thermal-structural problems, exact structural finite elements 
subjected to thermal loads are considered herein. 

3.3.1 Truss 

Typical thermal and structural models for truss elements are 
shown in Fig. 5. For a steady-state analysis, exact thermal finite 
elements for internal heat generation, surface convection and 
specified surface heating are presented in section 3.2. In this 
section Che exact element temperatures are used in Che development 
of truss elements for computations of displacements and thermal 
stresses. 

For a truss element subjected to a temperature change, thermal 
strain is Introduced in the stress-strain relation; 

^ " ^ref)l (3*50) 

where a„ is the axial stress, E Ut the modulus of elasticity, 
u is the axial displacement which varies with the axial coordinate 
X, a is the coefficient of thermal expansion, T(x) is the 
temperature, and Tj-gf is the reference temperature for zero stress. 
The rod equilibrium equation with an assumption of negligible body 
force is 



(3.51) 


which when combined with the stress-strain relation, Eq. (3.50), and 



+ Aa 
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Thermal and stress models of rod element 
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mulclpllod chrough by che truss cross-soctlonal area A, yields 
Che governing differential equation, 

EA^-oEA^ (3.52) 

Since the temperature T Is knovm from the thermal analysis, a 
general solution to the above differential equation can be obtained. 

An exact finite element can be formulated by applying the nodal 
displacements \x^ and U 2 as the boundary conditions at x ■ 0 
and X “ L, respectively. In this case, the exact element displace- 
ment variation Is 

X L 

u(x) - (a I T dx - a I f T dx) + (1 - J) u^^ 4= (|) u^ (3.53) 

0 0 


or In the matrix form 


u(x) - ^1^^^ ^^(x)] < 


u. 


r " l^ai (3.54a) 


where Nq(x) is the element nodeleas interpolation function,* 

N^, 1 » 1,2 are typical element interpolation functions, Uq is 
the nodeless parameter, and u^, 1 - 1,2 are the element nodal 
displacements. The element interpolation functions are 


X L 

Nq(x) - aj Tdx-a^j 
0 0 


T dx 


(3.54b) 
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N3^(x)-l--g N2 (x)-^ (3.54c) 

% 

where, for convenience, the nodeless parameter Uq is taken as unity 
in this case. Note that the element nodeless interpolation function, 
Nq(x), vanishes at nodes and depends on th<in integrals of element 
temperature variation obtained from the thermal analysis. 

To derive exact clement matrices, the method of weighted resid- 
uals is applied to the equilibrium equation (3.51). After performing 
an integration by parts and using the stress-strain relation, 

Eq. (3.50), element equations aind element matrices are obtained. 

These element equations are in the same form as those obtained from 
the variational principle described in section 2.3 and can be 
expressed as 

[Kg] {u} - {F^} (3.55) 

where [K ] is the structural element stiffness matrix, {u} is 
the vector of nodal displacements, and {F,j,} is the equivalent 
nodal thermal load vector. The elemei’t matrices are defined by (see 
Eqs. (2.27a) and (2.27d)) 

} dN dN 

^3..56a) 

0 

^ dN 

{F^} =» AEaj {-—} (T - Tj.^^) dx (3.56b) 

0 
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Using the exact displacGment interpolation functions, Eq, (3.54a), 
the element stiffness matrix above is a three by three matrix which 
contains coefficients , i,j -0,1,2. Since the governing 
differential equation, Eq. (3.52), can be cast in the self-adjoint 
form (see section 3.2), this element stiffness matrix is symmetric 
and Kq^, i “1,2 arc zero. Both the element stiffness matrix and 
the equivalent nodal thermal load vector can be evaluated in closed 
form as. 


CK3I 



{F^} 



(3.57a) 


(3.57b) 


where 


L 

“ AEa J (T - dx (3.57c) 

0 


Once exact nodal displacements are determined, exact displacement 
variation within an element can be computed from Eq. (3.53). Exact 
element stress can also be obtained b> substituting element displace- 
ment variation, Eq. (3.53), into the stress-strain relation, 

Eq. (3.50). In this case, the exact element stress in te.rms of 
nodal displacements is 
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(T - T^ef) dx 


(3. 58) 


Using the exact element temperature varlac:[•‘^ obtained from the 
thermal analysis (cases la-ld) , both element nodeless Interpolation 
functions Nq(x), Eq. (3.54b), and the equivalent nodal thermal load 
F^, Eq. (3.57c), can be evaluated In closed form as shown In Table 
3 and Appendix B, respectively. 


3.3.2 Hollow Cylinder 

For a hollow cylinder where the temperature T varies only In 
the radial direction (Fig. 6), the only non-zero displacement Is 
u(r) and all shearing stresses are zero. The radial stress 
and circumferential stress satisfy the equilibrium equation [22] 



(3.59) 


The stress-strain relations are 


\ » I + “(T - Tj-ef) (3.60a) 

% “ E + a(T - T^^f) (3.60b) 

®z “ E ^®z " ^'=^r + “ (T - Tref) (3.60c) 


where v is Poisson’s ratio; e^, Sg and are the radial, 

circumferential and longitudinal strain, respectively. 
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Table 3 

Truss Element Displacement Interpolation Functions, Ng(X)* 




Case 

Nq(X) 


Ifa^ 2 2 3 

— (X‘^-X) + — (-X + 3X - 2X'^) 


1(b) «{ 


T^-T^ cosh mL + TgCcosh mL-1) 


m SjLnli mli 


[(cosh mLX -1) 


T -T 

- X(cosh mL -D] + (sinh mLX - X sinh mL) } 


m 


1(c) 


a(T„-T.)L 2 9 

± ..— ^ — (X -X) + (-X + 3X 


2X^) 


1(d) 


^(Trt^T-)! 2 9 

-7 T ' - (X'^-X) + (-X + 3X^ 


2X^) 


*For all cases; N^_(X) = 1-X, N^(X) = X where X = x/L. 


OO 
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For the case of a thin hollow cylinder, the assumption of plane 
» 

stress (Og ■ 0) is used. Substituting the stress~strain relations, 
Eqs, (3.60a~b), into the equilibrium Eq. (3.59) and using the strain- 


displacement relations. 

e„ * 4^ and Ca » — (3.61) 

r dr • Or 

where u denotes the radial displacement, the governing differential 
equation for the case of plane stress is 

■ (l+v)af (3.62) 

A general solution to this’ differential equation is given by 

r 

u(r) “ (1 + v) j (T - T^ef) r dr + r + (3.63) 

0 . 

Since the radial temperature variation T is known from the 
thermal analysis (see section 3.2,2), the exact axisymmetric element 
displacement variation can be derived by applying the nodal displace- 
ments \i^ and U2 as the boundary cond'.tions at r « a and r » b, 
respectively. The exact element displacement variation is 

r 

u(r) = (1 + v) ^ j" (T - T^ef) r, dr 
a 

- (IH- V) I f <T- W 

^ (b2 - a2) ] 

a 



64 


•I* 


r a (b^ - r^) 
^ (b^ - a^) 



+ 


r ^ (r^ - a^) 

r (fa2 _ a2) 



(3.64) 


or in the matrix £om 


u(r) 


[NpCr) Nj_(r) 



LNsJ 


(3.65a) 


where NQ(r) is the element nodeless interpolation function; N^, 
i “1,2, are typical element interpolation functions, Uq is the 
nodeless parameter, and u^j^, i “1,2 are the element nodal displace- 
ments. The element interpolation functions are 


r 

N^(r) “ (1 + V) ^ I (T - Tref) r dr 

a 


- (1+ V) ^ 


(r2 - a2) 
(b2 - a2) 


b 

a 


(T - T^gf) r dr 


(3.65b) 


Ni(r) 


ri; (b^ - r^) , 
^ r (b2 - a2) ^ 


and 


»2(r) - [>- 




(3.65c) 


Like for the exact truss element, the nodeless parameter Uq is 
taken as unity, and the element nodeless parameter NQ(r) vanishes 
at nodes. Element matrices can be derived by following the same 
procedure described for the truss element. In this case, the element 
stiffness matrix and the equivalent nodal thermal load vector are 


b 
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[Kg] -J [Bg]'“^ [D] [Bg] r dr 

a 


(3.66a) 




b 

J [Bg]^ [D] {a} (T - Tygf) r dr 
a 


(3.66b) 


where [Bg] is the strain-displacement matrix obtained from 
Eq. (3.61), 




du 


"dNj_ 

dN2“ 


* 

^r 


dr 


"dx 

dr , 


X 

< 

^ m 4 


> m 



4 

> 

G ^ 


JU 





“2 



r 


r 

r 




[Bg] {u} 


[D] is the elasticity matrix (plane stress), 


(3.67a) 


[D] 



V 

1 


(3.67b) 


and {a} is the vector of coefficients of thermal expansion. 


(a) 



(3.67c) 


Using the exact element interpolation functions, Eq. (3.65), 
the element stiffness matrix and equivalent nodal load vectors 
corresponding to the heat transfer cases (see section 3.2.2) can be 
derived in closed form. Due to complexity of the element interpolation 
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functions, a computer-based s/mbollc manipulation language MACSVMA 
was used to perform the algebra and calculus required for these 
element matrix derivations. Results of these element matrices and 
exact element interpolation functions are shovm in Appendix B and 
Table 4, respectively. 

Once nodal displacements are computed, exact thermal stresses 
in both radial and circumferential directions can be determined. 
Using the stress-strain relations, Eq, (3.60), and the strain- 
displacement equations in the form of Eq. (3.67a), the element 
stresses can be written in terms of nodal displacements as 




r 

1 

> - [D] . 

r 

Q 

CD 

V. 


^ I"- 


[Bg] {u} - (1 + V) a (T - Tyef) \ 


(3.68) 


For the plane strain case (e^ = 0) , all equations formulated 
for the case of plane stress above may be used by replacing 
E.(l-v^) for E, v/(l-v) for v, and (l+v)a for a. In addi- 
tion, the longitudinal stress exists in this case and can be computed 
from the last equation of the stress-strain relations, Eq. (3.60c), 


°z “ " ^ref^ 


(3.69) 


3.4 Applications 

To demonstrate the capabilities of the exact thermal and 
structural finite elements developed in sections 3.2 -3,3, the finite 
element thermal analysis program TAP2 [23] and the finite structual 
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Table A 

Axisynanetric Element Displacement Interpolation Functions 



NflCr) 


-g+v) 

2rw 


UTi + ^ wTq) [r^inA » 


+ (T2 + WT.) [r^inCJS) . 

° " (b2-a2) ^ 

2b^ 


NjCr) 


a(b^-r^) 

r(b^-a^) 


N,(r) 


b(r^-a^) 

r(b^-a^) 
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analysis program STAP [24] are used. Elements discussed in this 
chapter were added to these programs. Conventional finite elements are 
also available in these programs, so comparisons between exact finite 

I 

elements and conventional finite elements could be made. 

3.4.1 Coffee Spoon Problem 

The exact rod element for conduction and convection described 
in section 3.2.1 is used for one-dimensional heat transfer in a 
coffee spoon [25], Fig. 7. The lower-half of the spoon submerged in 
coffee is convectively heated by the coffee at 339 K, and* the upper- 
half is convectively cooled by the atmosphere at a temperature of 
283 K. The ends of the one-dimensional spoon model are assumed to 
have negligible heat transfer. 

Three finite element models are used to represent the spoon: 

(1) two exact finite elements, (2) two linear conventional finite 
elements, and (3) ten linear conventional finite elements. Tempera- 
ture variations computed by these three finite element models are 
compared in Fig. 7. The figure shows that two conventional finite 
elements predict nodal temperatures with fair accuracy but are unable 
to provide details of the nonuniform temperature distribution includ- 
ing the zero temperature gradients at both ends of the spoon. The 
temperature variation obtained from ten conventional finite elements 
is in excellent agreement with the result froin two exact finite^ 
elements. It should be noted, however, that an approximate solution 
results from the use of conventional finite elements since the exact 
solution to the problem is in terms of hyperbolic functions, which 
were used in the exact element interpolation function (Appendix A, 

Case lb) . 
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Fig. 7. Conventional and exact finite element solutions 
for coffee spoon with conduction and convection. 
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3.4,2 Thermal StrooBOs In Hypersonic Wing 

A 136 member truss model of a hypersonic wing [26] » Fig. 3, was 
chosen to Illustrate the use of exact truss finite elements. The wing 
Is assumed to have varying convective heat along the leading edge, top 
and bottom surfaces and Is convectlvely cooled Internally. Tempera- 
tures along the wing root are specified. Two finite element thermal 
models are used to represent the wing truss. The first model 
consists of 136 exact condutlon-convectlon rod elements (see 
section 3.2,1) with one element per truss member. The second model 
Is Identical to the first model, but linear conventional finite 
elements are used, Fig. 9 shows a comparison of temperature distri- 
butions along the bottom members of the center rib of the wing truss. 
Results show that the exact finite element model provides a realistic 
temperature distribution which Iv characterized by higher temperatures 
near the center of each truss member and lower temperatures at the 
nodes. The conventional finite element model underestimates the 
actual temperatures and Is not capable of capturing the highly 
nonlinear temperature distribution along the rib. Therefore, further 
mesh refinement of the conventional finite element model Is needed 
if a realistic temperature distribution is to be predicted. 

For the structural analysis, both models employ the same 
discretization as in the thermal analysis. The structural boundary 
conditions consist of constraining the nodes along the wind root. 

Truss member temperatures obtained from the exact finite element 
thermal model are directly transferred to the exact finite element 
structural modeL for computations of displacements and thermal 
stresses (see section 3.3.1). Likewise, displacements and thermal 
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Fig. 8. Thermal structural truss model of a hypersonic wing. 
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Fig. 9. Comparison of temperature and stress distributions 
in wing truss, z«0. 
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stresses computed from the conventional finite element structural 
model are based upon linear member temperatures obtained from the 
conventional finite element thermal model. Comparison of the thermal 
stress distributions for the two analyses are made as shown in 
Fig. 9. The figure shows that conventional finite elements under- 
estimate member stresses with a relatively large error. This error 
is caused by the use of the inaccurate temperature distribution 
from the conventioniTii finite element thermal model. Comparative 
temperature and stress distributions of other wing sections (not 
shown) have similar trends. The results clearly demonstrate that 
Improved thermal-structural solutions can be obtained through the 


use of exact finite elements. 


Chapter 4 


MODIFICATION OF EXACT FINITE ELEMENT FORMULATION FOR 
ONE-DU-IENSIONAL LINEAR TRANSIENT PROBLEMS 

In the preceding chapter, exact thermal finite elements for one- 
dimensional steady-state heat transfer problems were presented. 
Steady-state element temperature interpolation functions were 
formulated ia closed form based upon solving ordinary differential 
equations. In transient analysis, exact element temperature inter- 
polation functions cannot be obtained in closed form since general 
solutions to typical transient problems are infinite series. However, 
by modifying the steady-state element temperature interpolation 
functions for the transient analysis, improved transient temperature 
solutions can be obtained as described in this chapter. 

In steady-state analysis, finite element temperature distribu- 
tions are a function of only the spatial coordinate, but for transient 
analysis, the element temperature distribution are a function of both 
space and time. For example, a one-dimensional transient heat conduc- 
tion is governed by the partial differential equation. 


RA 


3^T 

3x2 


pcA 


3T 

3t 


( 4 . 1 ) 


where k is the material thermal conductivity, p is the density, 
c is the specific heat, A is the conduction area, and T is the 
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temperature which varies with the spatial coordinate x and time t. 
A two-node linear conventional finite element may be us9d in the 
analysis where the element temperature variation T is expressed 
in the form (Fig. 10(a)) 


T(x,t> 




IjCt) 

I r 

T^(t) 




N2(x)J 


I 

T2(t) 


(4.2) 


where N^(x), i»l,2 are the element interpolation functions which 

are a function of the spatial coordinate x; L is the element length, 
and T^(t), 1*1,2 are the time-dependent nodal temperatures. 

With the heat equation shown in Eq. (4.1), the corresponding 
element equations and element matrices can be derived as described 
in section 2.3. Typical element equations have the form 


[C] {T} + [K] {T} * {Q} 


(4.3) 


where {T} and (T) denote vectors of nodal temperatures and the 

time rate of change of nodal temperatures, respectively. The matrix 

[k] and the vector {Q} represent the conductance matrix and the heat 

load vector, respectively, and have the same meaning as previously 

described for the steady-state analysis in the preceding chapter. 

The additional matrix [c] is called the capacitance matrix and defined 

* 

by (see Eq. (2.15a)) 


L 

[C] = j pcA {N^} dx 

0 


(4.4) 









For the linear interpolation functions shown in Eq. (4.2), the 
capacitance matrix can be evaluated as 
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[C] 


2 

1 


1 

2 


(4.5) 


This form of the capacitance matrix, Eq. (4.5), is called a consistent 
capacitance matrix because Its definition is consistent with the 
matrix formulation, Eq. (4.4). Quite often, the above capacitance 
matrix Is approximated by lumping the off-diagonal terms with the 
diagonal terms to give, 


[C] 


pcAL 


2 


1 


0 


(4.6) 


\ 

and is called a lumped capacitance matrix. It should be noted that 
degradation of the solution accuracy may result from the use of the 
lumped capacitance matrix compared with the consistent capacitance 
matrix. However, computational advantages (e.g. explicit time 
integration algorithms) may be achieved using the lumped capacitance 
matrix whereas the loss of solution accuracy may be insignificant 


[5]. 


4.1 The Nodeless Variable 

In the preceding chapter, exact finite elements for steady- 
state analysis are formulated based upon solving ordinary differential 
equations. Exact element temperature variations after imposing nodal 
temperatures as boundary conditions are written in the form, 
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T(x) - Nq(x) Tq + N^Cx) + N 2 (x) (4.7) 

where i *1,2 are element interpolation functions; T^, i ■1,2 

are unknown nodal temperatures, Nq(x) is the element nodeless 
interpolation function, and Tq is a known nodeless parameter. 

For transient thermal problems, it is not possible to formulate 
exact element interpolation functions in closed form because general 
solutions to typical transient problems are infinite series. However, 
since the transient response may approach exact steady-state solutions 
as time becomes large, the use of the exact steady-state element 
temperature variation in the form of Eq. (4.7) may provide better 
accuracy of solutions than those obtained from the linear conventional 
element, Eq. (4.2). 

To use the steady-state element temperature variation for 
transient analysis, Eq. (4.7) is written in the form, 

T(x,t) = NQ(k) Tq + Nj_(x) Tj_(t) + N 2 (x) T 2 (t) (4.8) 

where the unknown nodal temperatures T^ and T 2 become a function 
of time t. Since the nodeless parameter Tq is known and independent 
of time, the product of the nodeless interpolation function and the 
nodeless parameter, NqCx) Tq, retains the same shape throughout 
Che transient response. Characteristics of the element temperature 
variation expresr.ed by Eq, (4.8) during the response are illustrated 
in Fig. 10(b). 

Equation (4.8) may not be a good representation for a transient 
thermal responsv as will be shown by the following argument. As 
described in section 3.1.1, the nodeless parameter Tq is a scalar 
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quantity which contains a physical parameter associated with a given 
heat load. For example, the nodeless parameter for a slab subjected 
to a uniform internal heat generation rate is given by (see 

Table 2) 


where k is the material thermal conductivity and L is the 
element length. If a slab is modeled by an exact finite element, 
the element temperature variation is given by (see Appendix A, 
Case 2) 


L '' 


2k 


+ (1 - f) Ti + (f ) 


where and are the element nodal temperatures. If both 

surfaces of the slab have a specified temperature T in addition, 

s 

the above equation becomes 


T(x,t»0) 




+ (!-#) T, 


"s 


(4.9) 


For the case where the internal heat generation is raised instan- 
taneously from to Q^, the transient temperature variation 

within the slab should gradually increase and reaches the new steady- 
state temperature variation 


T (x , t-+“) 


i (1 _ i) ^ 
L ^ 2k 


+ (1 - f) T + (f) 

L S L 


(4.10) 
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where the new nodeless parameter is 


T 


0 



Since the nodal temperatures at both sides of the slab are fixed, 
the above two temperature variations, Eqs. (4.9 -4.10), suggest that 
the nodeless parameter Tq should vary with time so the element 
temperature can change gradually during the transient response. This 
argument leads to a modification of the element temperature variation 
employed in the steady-state analysis, Eq. (4.8), to the form 

T(x,t) - Nq(x) lQ(t) + Nj_(x) T^(t) + N 2 (x) T 2 (t) (4.11) 

where the nodeless parameter becomes an additional time-dependent 
element unknown and is called a nodeless variable. A typical element 
temperature variation with the nodeless variable TQ(t) is illus- 
trated in Fig. 10(c). 


4.2 Element Equations and Matrices 


in this section, element equations and matrices for both the 
nodeless parameter approach and nodeless variable approach are 
presented. In the nodeless parameter approach, the element tempera- 
ture variation shown in Eq. (4.8) can be written in the matrix form 



Tj_(t) ^ 

T2':t) 


T(x,t) = LNq(x) N^(x) N 2 (x)J 


< 


(4.12) 
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Using the definition of the element capacitance matrix shown in 
Eq. (4.4), the capacitance matrix for the given element interpolation 
functions can be derived. Element equations, Eq. (4.3), can then be 
written explicitly as 
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vrhere , i,j *0,1,2 are typical terms in the capacitance matrix; 

and Q^, i,j *0,1,2 are typical terms in the element stiffness 
matrix and the heat load vector previously described in the steady- 
state analysis. Since the nodeless parameter Tq is constant, its 
time-derivative 9TQ/3t is zero. Therefore, the first equation 
which involves the nodeless parameter is uncoupled from the nodal 
unknowns in the second and third equations. Hence, the above equa- 
tions reduce to 
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(4.14) 


Note that these equations contain two basic element nodal unknowns 
as for the linear conventional finite element. Once the nodal 
temperatures at a typical time are computed, element temperature 
variation can be obtained using Eq. (4.12). 

In the nodeless variable approach, the element temperature varia- 
tion shown in Eq. (4.11) can be written in the matrix form 
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To(t) 

where Tq(c) denotes the element nodeless variable which is a func- 
tion of time as the unknown nodal temperatures Tj^(t) and T2(t) . 
Since the element interpolation functions are identical to those 
used in the nodeless parameter approach, the element matrices are 
also identical. Element equations obtained using this approach have 
the form 

^0 
^1 
^2 

Because the nodeless variable is unknown, the equations are coupled 
through the capacitance matrix due to the presence of Tq. Thus 
typical element equations obtained from the nodeless variable approach 
contain three unknowns, l.e. one more unknown than the nodeless 
parameter apporach or the lineg^r conventional finite element. 

An advantage of the nodeless parameter and nodeless variable 
approaches is that both can provide an exact steady-state solution 
at the initial condition for the transient response. As time becomes 
large and new steady-state thermal equilibrium is reached, the exact 
temperature distribution may be predicted by both approaches. It 
should be noted that with the use of the nodeless variable approach 
the temperature variation within the nodeless variable element can 
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vary with Clmo even though the element nodal tcmperatucos are fixed. 
This feature is characterized by the term Nq(x) TQ(t) shown in 
Eq. (4.11) and is different from the linear conventional finite 
element where the clement temperature distribution is completely 
controlled by nodal temperatures. 

4,2.1 Rod Element 

The rod element with heat conduction combined with surface 
convection, Internal heat generation, or surface heating previously 
considered in Fig. 4, Cases la-d is extended for transient analysis. 
For each heat transfer case, the governing differencial equation for 
the temperature distribution T(x,t) can be derived using an energy 
balance on a small segment of the rod. These governing differential 
equations are: 


pcA II - IcA - 0 

3x^ 


(4.17a) 


pcA II - kA ^ + hpT - hpT, 
3x 


3^T 


(4.17b) 


peA II - kA ^ - QA 
3x 


(4.17c) 


pcA II - kA - qp 
3x'‘ 


(4.17d) 


where A is the element cross-section area, h is the convection 
coefficient, p is the cross-section perimeter, T^ is the surround- 
ing medium temperature, and q is the specified surface heating rate 


per unit area. As one example, conduatlon with Internal heat genera 
tlon where the exact stoady-state element temperature variation la 
(see Appendix A, Cane 1(b)) 
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I(=t,t) . f a - f) T;, + (1 - 2) Tj^ + (2) Ij 


LNq(x) Nj_(x) N 2 <x)J i Tjl > 


(4.18a) 


where Tq la the nodelesa parameter given by (see Table 2) 


^0 2k 


(4.18b) 


Using the exact element interpolation functions shown in Eq. (4.18a) 
above, the capacitance matrix is derived using Eq. (4.'f>.). The 
conductance matrix and heat load vector are derived using Eqs. (2.15b) 
and (2.16b), respectively. Therefore, the element equations are 
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> (4.19) 


In the nodeless parameter approach, the nodeless parameter Tq 


is constant and the above element equations reduce to 
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(4,20) 


wlch two unknown nodal temperatures Tj^(t) and l 2 (t), It should 
be noted that the element aqjjatlons obtained from using the nodeloss 
parameter approach shown in Eq, (4.20) above are identical to Chose 
obtained from Che linear conventional finite element for this heat 
transfer case. Thus, results of nodal temperatures during Che 
transient response are also identical. However, results of element 
temperatures are different due to the difference of their element 
interpolation functions, Eqs, (4.2) and (4.18b). As the transient 
response reaches the steady-state, the nodeless parameter approach 
provides exact solution for both nodal temperatures and element 
temperature variations where only exact nodal temperatures are 
obtained through the use of the linear conventional finite element. 

In the nodeless variable apporach, the element equations with 
two unknowns of nodal temperatures and an unknown of nodeless 
variable shown in Eq. (4.19) must be solved simultaneously. It can 
be seen from these equations that as the steady-state themal 
equilibrium is reached, the rate of change of nodal temperatures 
and the nodeless variable vanish. Then the first equation yields 


T 


0 



which is identical to the nodeless parameter shown in Eq. (4.18b). 
This means the nodeless variable varies during the transient response 
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and provides the value required for computation of the exact 
temperature variation when thermal equilibrium Is reached • 

For other heat transfer cases such as conduction with surface 
convection or surface heating, the same procedure Is applied. 

Element matrices corresponding to each heat transfer case in the form 
of Eq. (4.16) are given In Appendix C. Capabilities of the nodeless 
parameter and the nodeless variable finite elements for transient 
analysis are evaluated by comparisons with an exact transient 
conduction-convection solution and the linear conventional finite 
element In the first example at the end of the chapter. 


4.2.2 Axisymmetrlc Element 

Similar to the rod element, the axisymmetrlc element previously 
described in the steady-state heat transfer (Fig. 4, Case 3) is 
extended for the transient analysis. Radial heat conduction is 
combined with internal heat generation and specified surface heating 
or surface convection on the inner or outer cylinder surfaces are 
considered through the boundary conditiono. The governing differ- 
ential equations for the cases of pure conduction and conduction 
combined with internal heat generation are 


3t r 3r ^ 3r^ 


(4.21a) 


pc 


3T 

at 


k _a_ 
r 3r 


(r 


Q 


(4.21b) 


respectively, where r denotes the radial coordinate. 

Element equations can be derived by the method of weighted 
residuals applied to the governing differential equations. 

^ '<L 
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Typical oloraont aquations arc in the form of equation (A. 16), The 
element conductance matrix and heat load vector are identical to those 
obtained in the steady-state analysis shown in Eqs. (3,3db) and 
(3.38c), respectively. The element capacitance matrix associated with 
the rate of change of nodal topporaturea has the fonin 

b 

[C] - j pc {N^} [N^J r dr (4.32) 

a 

For the element interpolation functions associated with the heat 
transfer cases shown in Appendix A, the corresponding capacitance 
matrix can bo evaluated in closed form. Capacitance matrices in 
the form of Eq. (4,16) are given in Appendix C. 

4.3 Applications 

4,3,1 Transient Heat Conduction in a Rod with Surface Convection 
A rod with length L subjected to surface convection and 
specified end temperatures is sho^^n in Fig. 11(a). Initially the 
rod is convectively cooled by a surrounding temperature at 255 K 
and, at time t » 0”^, the s''t,rounding temperature is raised 
instantaneously to 589 K. Transient temperature distributions along 
the rod are computed using: (1) the exact solution [27], (2) two 

linear conventional finite elements, (3) two nodeless parameter 
finite elements, and (4) two nodelesa variable finite elements. In 
each finite element model, the element lengths are taken to be equal 
(L/2) with an unknown of nodal temperature at the center of the rod. 
Comparisons of the temperature variations obtained from these three 
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(a) Rod heated by surface convection. 


Fig. 11. Conventional and nodeless variable finite element 
solutions for a rod with surface convection. 
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flnitG olcmenc models and the exact solution £or t * 0, 0.01 and 

0.3 s, are made as shown in Fig. ll(b-d). 

At time t ■ 0 while the rod is in thermal equilibrium, two 
nodeless variable finite elements provide the exact steady-state 
temperature distribution. Two conventional finite elements are 
unable to provide details of the nonuniform temperature distribution 
due to the assumption of linear temperature distribution within the 
element. At time t ■ 0.01 s. (Fig. 11(c)) after the rod has been 
conveotively heated, the differences in the transient response 
predicted by three finite element models are shown clearly. Two 
linear conventional finite elements predict the unknown nodal 
temperature at the center of the rod with fair accuracy but the 
element temperature distributions are overestimated from the actual 
temperature distribution with a relatively high error. Two nodeless 
parameter finite elements yield the unknown nodal temperature with 
the same accuracy as of two linear conventional finite elements but 
predict extremely poor element temperature distributions. Two nodeless 
variable finite elements provide the best approximation of the 
unknown nodal temperature ^d.th excellent temperature distributions 
within the elements. As the rod temperatures approach a new steady- 
state solution at time t ■ 0.3 s. (Fig. 11(d)), two linear 
conventional finite elements yield a fair approximation of the 
unknown nodal temperature but crudely approximate the temperature 
distribution. Both the nodeless parameter finite elements and the 
nodeloss variable finite elements provide excellent prediction of 
the unknoTO nodal temperature and details of the nonuniform tempera- 


ture distribution. 
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(b) Comparative temperature distributions at t »0 s. 


Fig. 11. Conventional and nodeless variable finite element 
solutions for a rod with surface convection. 
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(c) Comparative temperature distributions at t ».01 s. 


11. Conventional and nodeless variable finite element 
solutions for a rod with surface convection. 
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(d) Comparative temperature distributions at t ».30 s. 


Fig# 11, Conventional and nodeless variable finite element 
solutions for a rod with surface convection. 
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Even chough the nodelcss paramecer approach employed in this 
problem yields excellent representation of the temperature distribu" 
tions at the beginning (t ■ Os,) and near the end (t ■ 3.0 s.) 
of the response, the approach is unable to provide reasonable element 
temperature distributions during the response. As shown in Fig. 11(c) 
at time t ■ 0.01 s., temperatures obtained from the nodeless 
paramecer finite elements are characterized by bumps within the 
elements. These unacceptable results are caused by using the steady- 
state element temperature distribution wi,th the constant nodeless 
parameter for Che transient analysis. Therefore, the nodeless 
parameter approach should not be employed for transient response 
predictions. Instead, the nodeless variable approach should be used 
since it gives acsuracy superior to the linear conventional finite 
elemenu throughout the response and predicts exact steady-state 
solutions . 

4.3.2 Transient Thermal Stresses in a Rod with Internal Heat 
Generation 

To further illustrate the use of the nodeiess variable approach 
for one-dimensional transient problems and demonstrate additional 

I 

benefits that can be achieved, an analysis of transient thermal 
stresses in a rod with internal heat generation is presented. 

A rod with constant, cross-sectional area A and length L 
encased between fixed walls is shown in Fig. 12(a). Both ends of 
the rod have the specified temperatures at 311 K and 333 K at 
X ■ 0 and x “ L, respectively. Initially, the rod is subjected 

3 

to a uniform internal heat generation rate Q » 353 kW/m and is 
in the thermal equilibrium. At time t ■ 0"^, 
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(a) ROD WITH INTERNAL HEAT GENERATION 


Fig. 12. Conventional and nodeless variable finite 
element solutions for a fixed end rod with 
internal heat generation. 


0 

533 K 


9S 


3 

generation rate increasea abruptly to 1073 kW/m and remains 

constant thcreaftor. The rod is modeled using: (1) 20 linear 

conventional finite elements, (2) two linear conventional finite 

elements, and (3) two nodeless variable finite elements. Comparative 

temperature distributions at time t •• 0, 0.1 and 1.0 hr. are 

♦ 

shown in Fig. 12(b). The figure shows that two nodeless variable 
finite elements have the same cc.{.*>ability in predicting transient 
temperatures as 20 linear conventional finite elements. Two linear 
conventional finite elements underestimate the temperature distribu- 
tions with relatively large error throughout the transient response. 

In the structural analysis, three structural finite element 
nodels with the same discretizations as for the thermal finite 
element models are employed. Element temperatures obtained from 
the thermal finite element model are transferred directly to the 
structural finite element model for computation of displacements 
and stresses. For the quasi-static analysis, the structural response 
are computed at times corresponding to the transient thermal solu- 
tions obtained previously. At each time, the equivalent nodal 
thermal forces are computed using Eq. (3.57) and the element nodal 
displacements are computed from Eq. (3,55). Once the element nodal 
displacements are obtained, element displacement distributions and 
element thermal stresses are computed from Eqs. (3.53) and (3.58), 
respectively. 

Displacement distributions obtained from the three structural 
finite element models are shown in Fig. 12(c). The figure shows 
that two linear conventional finite elements are inadequate to 
represent the details of nonuniform displacement distributions. 
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(b) COMPARATIVE TEMPERATURE DISTRIBUTIONS 


Fig. 12. Conventional and nodeless variable finite 
element solutions for a fixed end rod with 
internal heat generation. 
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(c) COMPARATIVE DISPLACEMENT DISTRIBUTIONS 


Fig. 12. Conventional and nodeless variable finite 
element solutions for a fixed end rod with 
internal heat generation. 
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Similar Co cho ehp.rmal analyoio, displacomonC dlscrlbutions obcained 
from cwo nodolosa variable flnlCQ olomenCs and 20 linear eonvonc:.onal 
finite olomonCs are in excolxont agroomenc throughout the transient 
response. Comparative thermal stresses obtained from these finite 
element models at the times mentioned above are given in Table 5. 
Thermal stresses computed fiom two nodelosr. variable finite elements 
and 20 linear conventional finite elements are equal since temperature 
variations of these two finite element models coincided. Two linear 
conventional finite elements underestimate the thermal stresses and 
the error Increases with time with a maximum of 10% at t ■ 1.0 hr. 

These two examples clearly demonstrate benefits of using the 
nodeless variable approach in one-dimensional transient thermal- 
structural problems. Further applications of the nodeless variable 
approach can be found in Ref. [28 1. The use of the nodeless variable 
for improving temperature solutions in the transient thermal 
analysis directly Improves accuracy of displacement and stress 
distributions in the structural analysis. The advantages of the 
nodeless variable approach for linear transient thermal-structural 
problems have been demonstrated in this chapter. The approach will 
be extended to nonlinear steady-state and transient thermal-structural 
analysis which includes radiation heat transfer in the next chapter. 


’ z' 
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Table 5 



Comparative Thermal Stresses for a Rod with 
Internal Heat Generation 

9 




Stress, HFa 


Time, t 
Hr. 

2 Conventional 
Elements 

20 Conv. Elements 
2 Nodeless Variable 
Elements 

% Diff. 


0 

-507 

-531 

4.5% 

1 

0.1 

-598 

-640 

6.5% 

1.0 

-652 

-724 

10.0% 


Chapter 5 


ONE-DIMENSIONAL THERMAL-STRUCTURAL FINITE ELEMENT ANALYSIS 
WITH RADIATION HEAT TRANSFER 

Due to their relatively low weight, high stiffness and ease of 

» 

fabrication, trusses have high potential for use In space structures 
for solar collectors, antenna and space stations i Thermal anal}^sls 
of these structures Includes conduction heat transfer combined with 
significant radiation heat transfer. Radiation heat transfer Intro- 
duces a strong nonlinearity In the energy equation being solved. 
Furthermore, a time dependent solution procedure Is required for the 
analysis due to the changing orientation of the structure during the 
orbit. 

In this chapter, finite element solution procedures for one- 
dimensional transient thermal analysis with radiation heat transfer 
are presented. Three finite element types are formulated: (1) an 

isothermal element, (2) a linear conventional element, and (3) a 
nodeless variable element. Accuracy and efficiency of the finite 
elements are evaluated using two thermal-structural examples at the 
end of the chapter. 

5.1. Solution Procedures for One»-dimensional Transient Thermal 
Analysis with Radiation Heat Transfer 

In this section, transient thermal analysis for a one- 
dimensional finite element with radiation heat transfer is presented. 


100 



101 


ThQ radiation surface Is assumed to be diffuse, gray and opaque 
which means the emitted radiation energy is uniformly distributed, 
independent of wave length and the material does not transmit 
radiation. For convenience all material thermal properties are 
assumed constant. 

For one-dimensional transient heat conduction in a rod with 
surface radiation, the governing differential equation for the 
temperature distribution T(x„t) can be derived using an energy 
balance on a small segment. With the assumptions mentioned above, 
the governing differential equation is 

2 

pcA 1^ - kA + e opg - a p q*^ (5.1) 

where p is the density, c is the specific heat, A is the rod 
cross-sectional area, k is the material thermal conductivity, 

0 is the Stefan-Boltzmann constant, e is the surface emissivity, 
a is the surface absorptivity, is the incident surface heating 

rate from distance directional sources per unit area, p^ and p^ 
are the cross-sectional perimeters for surface emitted energy and 
incident energy, respectively. 

Finite element equations corresponding to the governing differ- 
ential equation (5.1) can be derived using the method of weighted 
residuals as described in section 2.3. For this case, typical 
element equations have the form 

[C] {T} + [[Kj + [Kp]] {T} = {Qc> + {Q^} (5.2) 


where [c] is the element capacitance matrix; [K^] and [K^] 


are 


102 


the element conductance matrices corresponding to conduction and 
radiation, respectively, {Q } Is the element vector of conduction 

w 

heat flux across element boundary, and {Q^} Is the element heat 
load vector due to Incident radiation. These matrices are expressed 
In the form of Integrals over the element length L as follows: 


L 

[C] * j* PcA {Nrj,} [NijJ dx 

0 
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[K^] {T} « J 6 oPg dx 
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L 

{Qj,} - J ap^ {N^} dx 

0 


(5.3a) 


(5.3b) 


(5.3c) 


(5.4a) 


(5.4b) 


where [n^J denotes the element temperature Interpolation functions. 

As shown In equation (5.3c), the conductance radiation matrix 
contains the element temperature within the Integral. The element 
equations, Eq. (5.2), thus constitute a nonlinear set of equations. 
Since the time rate of change of the temperature vector {t} also 
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appears In the element equations, a transient nonli. .ar solution 
procedure Is required for the analysis. 

Typical techniques for transient nonlinear solutions combine 
a linear transient solution method and a steady-state nonlinear 
solution method. The solution technique here uses a tlme-marchlng 
scheme where temperatures are computed at the middle of the time 
step, the Crank-Nlcolson algorithm. At each time step, Newton-Raphson 
Iteration Is used to correct for nonllnearltles . Further details of 
these methods Including other solution algorithms can be found In 
the finite element text, Ref. [l5]. 

Starting from the element equations, Eq. (5.2), the tlme- 
marchlng scheme Is first applied by approximating the time rate of 
change of nodal temperatures as 

(t) -jt - (T>„) (5.5) 

where At Is the time interval between the time step n and n+1 

such that t , 1 » t + At: {T} and {T} ,, are the vectors of 

n+1 n n n+1 

nodal temperatures at the time step n and n+1, respectively. 

Since the Crank-Nlcolson algorithm computes temperature solutions at 
the middle of time steps, nodal temperatures at the middle of the 
step are approximated by 

{T} - I ({T}^ + (5.6) 

vhere {T} denotes the vector of nodal temperatures at the middle 
of the step. From this equation, the vector of nodal temperatures 
at the step n+1 is 
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- 2CT} - (T}j^ (3.7) 

By combining Eqs. (5.5) and (5.7), the time rate of change of the 
nodal temperature vector shown in Eq. (5.5) can be expressed in 
terms of nodal temperatures at the middle of the step and step n 
as 


(T) - ({T} - {T)n) (5.8) 

Substituting Eq. (5.8) into Eq. (5.2), the element equations become, 
[C] + [K^] + [K^]] {T> - {Q^} + {q,.} [C] {T}^ (5.9) 

In Eq. (5.9), the vector of nodal temperatures (T}^ that 
appears on the right-hand side is known from the previous step. 

Since the unknown nodal temperatures contained in the vector {T} 
are computed at the middle of time step, the heat load vectors must 
be evaluated at the same time. Once the unknown nodal temperature 
vector (T) is obtained, the nodal temperature vector at 

the step n+1 can be computed from Eq. (5.7). 

The element equations obtained by applying the Crank-Nicolson 
algorithm shown in Eq. (5.9) are in the form of nonlinear algebraic 
equations 

CK(T)] {T} - {q} . (5.10) 

where 

[K(T)3 {T} - [-^ [G] + [KJ + [Kj.]] {T} 


(5.11a) 
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and 

{Q} - (Qc) + {Q^} [C] {T}^ (5.11b) 

For any temperature vector (x) that Is not an exact solution to 
equations shown In Eq. (5.10) above, an unbalanced nodal heat loads 
exist which can be written In the form of a vector {’!') as 

{t|/} « [K(T)] {T} - {Q} (5.12a) 

or In tensor notations, 


r 

- Z 

j-1 


K 


Ij 


Tj -Q, 


(5.12b) 


To develop the Newton-Raphson method a Taylor aeries exp an sion, with 
the first order-derivative accuracy Is written as 


'i»i({T}®) + ■( )“ ATj 


(5.13) 


A set of algebraic equations Is obtained In the form 

[J]“ - {R}“ (5.14) 

where the superscript m denotes the m*^^ Iteration. The matrices 
[J]“ and {R}“ are the Jacobian matrix and the residual load vector, 
respectively, defined by 

^11 'Bj 


(5.15a) 
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\ ’^1 (5.15b) 

At each iteration, the vector of nodal temperature increments 
is computed using Eq. (5.14) and a new temperature vector 
is obtained from 

« {T}“ + (5.16) 

The iteration process is terminated when a convergence criteria 
(such as the maximum nodal temperature increment is leas than a 
specified value) is met. For steady-state analysis, the equations 
shown in Eq. (5.2) do not contain the time rate of change of nodal 
temperatures, and only the Newton-Raphson iteration is required. 

5.2 Element Formulations 

In this section, three one-dimensional finite elements with 
Surface radiation are formulated. Crank-Nicolson and Newton-Raphson 
methods described in the preceding section are employed for the 
transient and nonlinear solutions, respectively. 

5,2.3. Isothermal Element 

The isothermal element is a simple finite element suitable for 
problems with negligible conduction heat transfer. A uniform 
temperature variation is assumed along the element that varies only 
with time (Fig. 13(a)). This element is different from the finite 
elements mentioned in the previous chapters since element temperature 


is the only unknown for the isothermal element, Since the element 





(a) Isothermal element 




Fig. 13. One-dimensional eleaent interpolation 
functions for nonlinear transient 
analysis with radiation. 


neglects heat conduction, the governing differential equation shown 
In Eq. (5.1) becomes 
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dT 4 

pcA^+ eopg T - apq p^. 


(5.17) 


where T denotes the element temperature which Is a function only 
of time t . 

The CranU-Nlcolson algorithm lo applied to the above differential 
equation by first writing the rate of change of the element tempera- 
ture In the form of Eq. (5.8), 


f ■ 


where T denotes the element temperature at the middle of the step. 
Substituting this equation Into Eq. (5.17) yields a nonlinear 
algebraic equation In the form 

(~ pcA + e opg T^) T - a p^ q^ + ~ pcA (5.18) 

After the element temperature T at the middle of the step shown in 
the above equation is obtained, the element temperature at the end 
of the step is computed from 


T 


n+1 


2T - T 

n 


C5.19) 


Next, the nonlinear algebraic equation shown in Eq. (5.18) is 
solved by applying the Newton-Raphson method. In this case, the 
unbalanced element heat load is given by (see Eq. (5.12)), 





109 

tp " ( ^‘ + e apg T^) T - a Pq ~ s^cA 7^ (S. 20 ) 

Using Taylor series approximation, an algebraic equation is obtained 
in the form 


j“ - R® (5.21) 

where J® and R® are the Jacobian and the residual load at the 
iteration defined by 

" It " St ^ (5.22a) 

R® - ~ PcA (T^ - T) - so Pg T^ 

(5.22b) 

+ a Pq Pr 

At each iteration, the element temperature increment AT is computed 
from Eq. (5.21) and a new element temperature is obtained from 



After the convergence criterion is met, the element temperature 
shown in Eq. (5.23) is used in Eq. (5.19) to compute the element 
temperature at the end of the step, The use of the isothermal element 
does not require a set of simultaneous equations due to the assump- 
tion of negligible heat conduction as previously mentioned. The 

< 

transient response of each element, therefore, can be computed 
separately. 
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The loothormal elemone la uaeful for modoling Crusa memboro 
whoro boat: conduction lo ncgllglblo in coroparloon with the Incident 
heating and omitted radiation. Applications of the iaothermal 
olcraont for transient onalysio of truss-type structures with surface 
radiation can be found in Refs. 129»30]. 


5.2,2 Conventional Element 

For the conventional element, a linear temperature variation 
is assumed between the two element nodes (Fig. 13(b)), 


T(x,t) 




[ T^Ct) 


Ijd:) 




(3.24) 


where the unknown nodal temperatures T^^Ct) and T 2 (t) are a 
function of time t. With the conduction-radiation differential 
equation shown in Eq. (5,1), element equations esn be derived as 
shown in Eq. (5.2). The vector of unbalance nodal heat loads show,i 
in Eq. (5.12a) is written explicitly in the form 


Cl'} [C](T} + - {Q^} - {q^} 


J. 

At 


[C]{T}^ 


or 


{'/'} - [g]{(t} - (T) } + [K Ht} + 

iiC net 

- {Q,} - {Q^} 


(5.25) 
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Foe sonvonionco nodal hoac load veetora coerooponding to each tom 
on tho right-hand oldo of tbo above cquotlon are introduced to yield 

W " <5.26) 

The Jacobian matrices and the residual heat load vector can now be 
formulated by using the definitions (see Eq. (5.15)) 

■'tj ■ ^ *1 ■ - 

For example, the first term on the right-hand side of Eq. (5.26) is 
the heat load vector associated with the capacitance matrix, 


'V tcHW = {Ty 


L 

- J pcA {N^} [N^J dx {{T} - {T}^} 
0 


With the Linear element interpolation functions shown in Eq. (5.24), 
this term can be evaluated in closed fom as 



Using the definition of Jacobian Jj^j ■ 3i|/^/3T^, i,j ■ 1,2, the 
corresponding Jacobian matrix is 
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(3.27) 


Similarly* cho Jacohlan maCrix aasociatod with the conductance conduc- 
tion matrix la obtained In the form 






(5.28) 


The third term on the right-hand side of Eq. (5.26) Is the heat load 
vector associated with the radiation matrix, 


L 

- [Kj] (T) •» f eo Pg {N-j} dx 
■^0 
or 

L 

® ^ Pg dx 

Therefore, the corresponding Jacobian Is 


■"ij ■ 3i:' M ” "i "j 
^ 0 

I 

or 

* 

L 

®°^Ps dx (5.29) 

0 
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With the linear element Interpolation functions shown In Eq. (5.24), 
this Jacobian matrix can be evaluated In closed form as 




15 


eaPgL 


10tJ+6tJt2 +3T^t| +t| 

2tJ +3tJT2 +3Tj^T^ +2T| 


2tJ +3tJT2 +3Tj_t| +2t| 

tJ +3tJT2 +6Tj_t2 +10T| 

(5.30) 


It can be seen that the Jacobian matrix associated with the 
radiation conductance matrix Is strongly nonlinear since the unknown 
nodal temperatures contribute to all terms In the matrix. The matrix 
Is sometimes [31] approximated by lumping these terms together 
similar to the lumped capacitance matrix given In Eq. (4.6). The 
lumped Jacobian matrix results in a much simpler form with zero off- 
diagonal terms, 



2 eopg L 




(5.31) 


From Eqs. (5.15) and (5.26), the total residual load vector 
is 


{R} - - 


{<PC^} - + {>p 



(5.32) 


For example, the residual load vector associated with the radiation 


conductance matrix is 
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{■Rk^^ “ tT} 

L 

" - \ e a Ps {N^} dx 

i) 

Using the linear element interpolation functions shown in Eq. (5.24), 
this residual load vector can be evaluated in closed form, 

5tJ + 4tJt2 + 

!■ j 

+ 2tJt2 + 3T^t| + 4Tj_t| + 5tJ 

After all Jacobian matrices and residual heat load vectors are 
computed, a final set of algebraic equations is obtained in the form 

[J]“ {AT}®^^ - {R}“ (5.34a) 

where the superscript m denotes the mth iteration and, 

[J] - [Jcg] + [JKg] + (5.34b) 

m - {RCc> + ^%r^ t^c> + ^^r^ (5.3^c) 

The solution of the temperature vector at successive times proceeds 
as previously discussed for the isothermal element. 

As shown in Eq. (5.34a), the transient and nonlinear solution 
procedures lead to a set of algebraic equations. The Jacobian 
matrices and the residual heat load vectors shown in Eqs. (5.34b-c) 
are thus necessary for the analysis.. With the linear element inter- 
polation functions shown in Eq. (5.24), these matrices can be 
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evaluated In closed form and are giv >n as computer subroutines In 
Appendix D. 

5.2,3 Modeless Variable Element 

In the preceding chapter, the nodeless variable approach was 
introduced for improvement of temperature solutions in one-dimensional 
linear transient analysis. The basic idea of the approach is the use 
of the steady-state element temperature intei'polation function derived 
from the solution of a given ordinary differential equation. An 
element nodeless variable is employed so that exact steady-state 
solutions are obtained at the beginning and at the end the 
transient with realistic temperature distributions prediced throughout 
the response. For one-dimensional conduction-radiation heat transfer, 
it is not possible to obtain a closed form solution to the governing 
differential equation, Eq. (5.1). However, the nodeless variable 
approach is still useful for the analysis to provide improved 
temperature solutions for the thermal element while maintaining the 
same discretization as the two node structural element. The element 
temperature distribution with a nodeless variable is written in the 
form, 

T(x,t) - LNq(x) N^(x) N 2 (x)J Tj_(t) y {T} (5,35) 

T2(t) 

^ - 

where Nq(x) is the nodeless variable interpolation function, 

N^(x), i =1,2 are typical element interpolation functions; TQ(t) 
is the nodeless variable, and T^(t), i =1,2 are the nodal tempera- 


tures'. 


U6 


As mentioned earlier, the nodeless variable Interpolation £unc** 
tlon Nq must vanish at nodes In order to preserve continuity of 
temperaturf/ between elements. There are wide choices for selecting 
the nodeless variable Interpolation function to meet »'hls requirement. 
The simplest function Is In the fom of polynomials with one order 
higher than the linear element Interpolation functions used in the 
conventional finite element, 


- ! <1 - f> 

(5.36a) 

Hj(x) - 1 - ^ 

(5.36b) 

Nj(x) - f 

(5.36c) 


With these element interpolation functions, the element temperature 
distribution, Eq. (5.35), results in a parabolic distribution over 
the element length as illustrated in Fig. 13(c). 

The use of the nodeless variable elero,ent for transient heat 
conduction with surface radiation follows the same procedure 
described for the linear conventional element. Typical element 
equations derived from the method of weighted residuals shown in 
Eq. (5.2) contain three unknowns. These element unknowns are the 
nodeless variable T^ and two nodal temperatures T^^ and 
For transient solutions, the Crank-Nicolson algorithm is applied, 
and a set of nonlinear algebraic equations is obtained. Next the 
Newton-Raphson method is used and a new form of simultaneous 
Eq. (5.34a) is obtained where the Jacobian matrices and the residual 
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load vectors are defined by Eqs, (3.54b"c), respectively. For example, 
the Jacobian matrix contributed from the radiation conductance matrix 
has the form 


L 

[J^ ] - 4 f e a Pg {N^} dx 
^ '*0 


Using the element nodeless variable temperature distribution and 
their element interpolation functions shown in Eqs, (5,35) and (5.36), 
respectively, this Jacobian matrix is written explicitly as, 

[Nq n^J dx 


[J 1 - 4 f EC (N„Tg + + BjTj)' 


N, 




Due to the complexity of the Jacobian matrix as shown above and 
other matrices that appear in Eq. (5.34), the computer-based symbolic 
manipulation language MACSYHA [32] was used to evaluate the matrices 
in closed fom. Results of the Jacobian matrices and the residual 
load vectors are provided in the form of computer subroutines in 
Appendix D. 

After the Jacobian matrices and the residual load vectors are 
computed, typical element equations shown in Eq. (5.34a) can be 
written in the form, 



m 

/• \ 

m+l 


m 

CM 

O 

o 

o 

o 

1 



^0 


10 “^ll “^12 

•< 


^ S3 ^ 

^1 

> (5.37) 

20 '^21 *^22 


ro 


^2 



V. J 


(5.37) 
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These element equations contain unknowns In the Increment^^f of the 
nodelss variable and the nodal temperatures, l.e. one more unknown 
than those obtained from the linear conventional element , Once these 
unknowns are obtained, new values of the nodal temperatures and the 
nodeless variable are computed using Eq. (5.16). After the Iteration 
process Is terminated, the nodal temperatures and the nodeless 
variable at the end of time step are computed from Eq. (5.7). 

Finally, the temperature distribution within the element Is computed 
by using the element nodeless variable Interpolation functions shown 

t 

In Eq. (5.35). 

It should be noted that the nodeless variable Interpolation 
functions, Eq. (5.36), Introduced In this section are applicable 
when other heat transfer modes (such as surface convection) are 
Included In the analysis. The element temperature distribution In 
the parabolic form can provide a more realistic .temperature dlstribu 
tlon than the linear conventional element. This type of the nodeless 
variable interpolation functions suggests that the nodeless variable 
approach can be generalized to other finite element types. To 
investigate this possibility, a two-dimensional nodeless variable 
thermal element is developed in the next chapter. 

5.3 Applications 

The effectiveness of the nodeless variable finite element 
described Jn this chapter is demonstrated for two examples of conduc- 
tion and radiation heat transfer. The linear conventional finite 
element described in section 5.2.2 is used in these two examples for 
comparison of solution accuracy. Temperatures computed from the 
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nodeless variable and linear conventional finite elements are used 
in the structural analysis for computation of displacements and 
thermal stresses. 

5.3.1 Thennal Stress in a Rod with Surface Radiation 

A rod with constant cross-section area A and length L 
encased between fixed walls is shown in Fig. 14(a). The rod has 
specified end temperatures at 311 K and 533 K at x ■ 0 and 
X ■ L, respectively, and is cooled along the surface by radiation 
to zero medium temperature. The rod is modeled using (1) 20 conven- 
tional elements with consistent Jacobian matrices (see Eq. 5.30)), 

(2) txro conventional elements with consistent Jacobian matrices, 

(3) two conventional elements with lumped Jacobian matrices (see 
Eq. (5.31)), and (4) two nodeless variable elements. The terms 
consistent and lumped refer to the formulation of the Jacobian matrix 
contributed by the radiation conductance matrix described in section 
5.2.2. 

Temperature distributions computed from these four finite 
element models are compared as shorm in Fig. 14(b). The figure shows 
that two nodeless variable elements have the same capability in 
predicting the unknown nodal temperature (at x/L ■ 0,5) and element 
temperature distributions as 20 conventional finite elements. Two 
conventional finite elements with consistent formulation underesti- 
mate the unknown nodal temperature and crudely approx 'mate temperature 
distribution. Two conventional finite elements with lumped formula- 
tion overestimate both the unknown nodal temperature and element 
temperature distributions with relatively large error. 
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RADIATING HEAT 



(a) ROD WITH SURFACE RADIATION 


Fig. 14. Conventional and nodeless variable finite 
element solutions for a fixed end rod 
radiating to space. 




X/L 

(b) COMPARATIVE TEMPERATURE DISTRIBUTIONS’ 


Fig. 14. Conventional nodeless variable finite element 
solutions for a fixed rod radiating to space. 
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In the atructural analysla, four structural elements with the 
same discretizations as for the thermal models are employed. Element 
temperatures obtained from the thermal model are transferred directly 
to the structural finite element model for computation of displace- 
ments and stresses. The conventional structural finite elements 

I 

employ linear element displacement distributions as used in typical 
finite element programs. The structural finite element for the 
nodeless variable thermal element uses the exact displacement 
distribution, Eq. (3.53), derived based upon the parabolic element 
temperature distribution shown in Eq. (5.35). Displacement distribu- 
tions obtained from these structural finite clement models are 
compared as shown in Fig. 14(c). The figure shows that two conven- 
tional finite elements are inadequate to represent the nonuniform of 
displacement distribution. In addition, two conventional finite 
elements with consistent and lumped formulations overestimate the 
thermal stress (not shown) by 12 and 23 percent, respectively. 
Displacement distributions obtained from two nodeless variable 
finite elements and 20 conventional finite elements are in excellent 
agreement where the difference in the thermal stresses is negligible 
(less than 0.05 percent). 

5,3.2 Thermal Analysis and Structural Response of a Space 
T russ Module 

A three member “orbiting truss module shown in Fig. 15 (a} is 
used to demonstrate the efficiency of the nodeless variable finite 
element. A typical truss member receives incident heating which is 
a combination of; (1) solar heating, (2) earth emitting heating, 
and (3) earth reflected solar heating. With the open-truss type 
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X/L 



(c) COMPARATIVE DISPLACEMENT DISTRIBUTIONS 


Fig. 14. Conventional and nodeless variable finite 
element solutions for a fixed end rod 
radiating to space. 
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otruecurc aa ahown in eho I’iguco, mombor Co mcmbor radlaclon oxebiingv^s 
are rclaelvoly small and arc neglccCed* A geoaynchronous orbit 
(period of 24 hr.) is employed where solar hoacing is large compared 
Co Che earch emlCCed heacing. During che orblc, incidenc heating 
normal Co a typical truss member varies continuously due to the 
changing orientation of the member. As the orbiting truss enters 
and leaves the earth's shadow, the incident heating changes rapidly. 
Member temperatures and structural deformations thus depend strongly 
on the time-dependent Incident heating. 

* 

To demonstrate the use of Che conventional and the nodeless 
variable finite elements formulated in the preceding section, the 
truss module with properties of aluminum Is considered. Four finite 
element models are employed where each truss member Is represented 
by: (1) 10 conventional elements with consistent formulation, (2) one 

conventional element with consistent formulation, (3) one conventional 
element with lumped formulation, and (4) one nodeless variable 
element. Temperature distributions computed from these four finite 
element models at a typical orbital position are shown In Fig. l5(b). 
The figure shows that the nodeless variable finite element model 
provides excellent pxrediction of the nodal temperatures and very 
good element temperature distributions compared to the refined 
conventional finite element model. The conventional finite elements 
with consistent formulation tend to average the nonuniform tempera- 
ture distributions and thus cannot provide accurate nodal temperatures. 
The conventional finite elements with lumped formulation predict 
nodal temperatures very well but yield large errors for member 
interior temperatures. Comparative temperature distributions of 
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finite element modela at other orbital positions aliio show a similar 
trend; the nodeless variable finite elements predict nodal tempera- 
tures and member temperature distributions very accurately compared 
to the refined conventional finite elements# 

Temperature obtained from the four thermal finite element models 
for a complete orbit are transferred to the structural finite element 
models for computation of displacement histories. The quasi-static 
analysis described in section 2.4 is employed for the computation 
of the unknown nodal displacements. Fig. 15(c) shows a comparison 
of typical member elongation histories computed from the finite 
element models during the orbit. Since the temperature distributions 
obtained from the nodeless variable finite element model and the 
refined conventional finite element model are in very good agreement, 
member elongation histories predicted by these two finite element 
models almost coincide (maxiiinim difference of 1 percent) . Conven- 
tional finite element models with consistent and lumped formulations 
yield errors for member elongation up to 29 and 44 percent, 
respectively. Such large errors result from the incapability of the 
conventional finite element to provide realistic member temperature 
distributions. Since the conventional finite element with consistent 
formulation trends to average the member temperature as previously 
mentioned, accuracy of the member elongation history obtained is thus 
higher than the conventional finite element with lumped formulation. 

These two examples clearly demonstrate the benefits of using 
the nodeless variable finite elements in one-dimensional radiation- 
conduction problems that are characterized by nonuniform temperature 
distributions. The elements predict member temperatures accurately 
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and are compaClble with two-node structural elements to permit an 

t 

integrated thermal-structural analysis. Additional applications of 
the nodeless variable finite element for one-dimensional thermal 
problems with conduction and radiation heat transfer appear in 
Ref. [33]. 


Chapter 6 


TWO-DK-IENSIONAL NODELESS VARIABLE FINITE ELEilENTS 

In the two preceding chapters the nodeless variable approach 
was applied to one-dimensional linear thermal- structural analysis 
and to nonlinear radiation heat transfer. The unique feature of the 
approach is the use of an additional nodeless variable for a thermal 
finite element. Improvement of solution accuracy is achieved while 
the same discretization is employed for both thermal and structural 
finite element models. 

In this chapter the approach is extended for development of 
two-dimensional nodeless variable finite elements. Restrictions for 
developing these finite elements are first discussed. Two nodeless 
variable finite elements and their interpolation functions are 
presented. Then the use of the nodeless finite elements for linear 
thermal-structural analysis is described. Efficiency of the nodeless 
variable finite elements is evaluated by comparison with the 
conventional bilinear four-node finite element and exact solutions 
in examples at the end of the chapter. 

For simplicity in understanding characteristics of the two- 
dimensional nodeless variable finite elements, a brief description 
of a conventional bilinear four-node thermal finite element is 
first given. The element temperature distribution for a bilinear 
four-node element is expressed in the form, 
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- LN-jJ tT} (6.1) 

where i “1,4 are the element interpolation functions which 

are a function of spatial coordinates in two-dimensions, and 
i -1,4 are the time dependent nodal temperatures. 

Fig. 16(a) shows a conventional four-node element with a general 
quadrilateral shape. As described in section 2.3, typical finite 
element matrices are in the form of integrals over the element 
volume or along the element boundary. Such element matrices for 
a quadrilateral shape are difficult to evaluate. To simplify the 
element integrations, the quadrilateral element in the Cartesian 
coordinate system (x,y) is transformed to a natural coordinate 
system (5,n) as shown in Fig. 16(b). The two coordinate systems 
are related by 

4 

X = T N,(£;,n) X. » InJ {x} (6.2a) 

i=l 

4 

y = E N. (5,n) y. = [nJ {y} (6.2b) 

i=l ^ 

where i »1,4 are the element shape functions defined by, 

Ni = -^(1-5) (1-n) ^2 “ 

N3 = jil+0 (l+n) = -|(l-5) (l+n) 


T-Uj_ N2 N3 N^J \ 


(6.3a) 
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(a) Global coordinates 


n 



(b) Natural coordinates 


Fig. 16. Four node isoparametric finite element in 
global and natural coordinates. 


133 


or In compact form, 

(1 + 55^) (1 + i -1,4 (6.3b) 

where 5^ and n^, 1-1,4 arc che nodal coordinates In the natural 

coordinate system. For example, ^2 " ^2 " 

etc. 

When the shape functions shown In Eq. (6.3) are used as 

the element temperature Interpolation functions In Eq. (6.1), this 
conventional element Is called an Isoparametric quadrilateral element 
because the some Interpolation functions are used to Interpolate 
temperature and spatial coordinates. 

Note that an element temperature Interpolation function shown 
In Eq. (6.3) has a value of unity at the node to which It pertains 
and a value of sero at the other nodes. Along the element edge 
(5 - ±1, n - ±1), these element Interpolation functions are linear. 
Therefore, the temperature distribution along a typical element 
edge varies linearly where the magnitude depends on the temperatures 
of the two corner nodes located at that edge. When elements are 
connected, the conventional quadrilateral element preserves 
continuity of temperature along the element interfaces. The conti- 
nuity of the element Interface temperatures is a basic requirement 
to assure convergence of the temperature solution as element size 
decreases. This continuity requirement must be met when a new 
thermal finite element is constructed. Further details of require- 
ments for a typical finite element to meet convergence criteria can 
be found in Ref. [15]. 
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6.1 Two-Dimensional Modeless Variable Thermal 
Finite Elements 

In several thermal-structural applications, a more detailed 

finite element thermal model is required than the finite element 

structural model. To maintain the same discretization for thermal 

and structural models, new thermal finite elements are required. In 

this section, two type of two-dimensional nodeless variable thermal 

finite elements for Improved temperature solutions are presented. 

These elements predict more realistic temperature distributions than 

♦ 

the conventional finite element previously described. The basic 
objectives for developing the new finite elements are: ( 1 ) the 

elements should provide a nonlinear temperature distribution but 
maintain four element nodes to be congruent with the four node 
structural element, and ( 2 ) compatibility of temperature along 
element interfaces must be preserved. The nodeless variable concept 
previously described for one-dimensional element is extended to 
two-dimensions to meet these objectives. 

6.1.1 "Bubble" Modeless Variable 

One approach. [34] for constructing nodeless variable finite 
elements is to add a "bubble" function which vanishes along the 
element boundaries. The element temperature distribution is written 
in the form, 

T - Ln^ N 3 nJ ^ 


Vo 


(6.4) 



135 


where Nq la the nodeloaa interpolation (bubble) function defined 
by 

Ng ■ (1 - C^)(l - nh (6.5) 

and Tq is the nodeless variable. 

Along the element boundary (5 •• ±1, o ■ 41), the bubble 
function, Eq. (6,5), is identically zero. Therefore, the element 
boundary temperature reduces to a linear variation as for the 
conventional finite element and continuity of temperature along 
element interfaces is preserved. Within the bubble nodeless variable 
element, the temperature distribution is a combination of the 
conventional element temperature distribution and a bubble function 
where its magnitude is measured by the nodeless varalble Tq. The 
combination thus permits a quadratic temperature distribution over 
the element. 

It should be noted Chat even though the bubble nodeless 
variable finite element can provide a quadratic temperature distribu- 
tion within the element, the temperature along the element boundary 
is linear. To achieve further improvement of the temperature solu- 
tion, the temperature distribution should vary nonlinearly along the 
element boundary. With the idea of the bubble function, a nodeless 
variable finite element with this behavior can be constructed. This 
type of nodeless variable finite element is presented in the next 


section. 
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6.1,2 Boundary NodolesB Varlablo 

Xn order to establish a nonlinear temperature distribution along 
the four element edges as well as within the element interiori the 
following four nodeless variable interpolation functions (see Fig, 17) 
are employed 


1 (1 - - n) 

(6.6a) 

1 (1 + oa - nS 

(6,6b) 

1 (1 - 5^)(l + n) 

(6.6c) 

4 (1 - ?) (1 - n^) 

(6 , 6d) 


where each interpolation function varies quadratlcally along one 

edge and vanishes on the other edges. For example, the nodeless 

2 

variable interpolation function varies as 1-5 along the 

edge n ■ -1 and is identically zero on the other three edges. 

As mentioned earlier, continuity of the temperature along the 
element interfaces must be assured for convergence of the solution. 
This restriction can be met by providing a nodeless variable for 
each element edge. With a nodeless variable for each element edge, 
element interpolation functions for a quadrilateral element can be 
written in the form, 



"1 



"5 

T - N2 N3 Nj ^ 

"^2 

■ + L »5 

»7 “aJ ' 

^6 

"7 


^4 

^ > 



^8 


( 6 . 7 ) 


OF POOR QUALITY 
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"7 ■ <1+") “S ■ 


Fig. 17. Nodeless variable interpolation functions for 
two-dimensional quadrilateral finite element. 
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whore T^i 1 "IjA and i "5, 8 are the nodal tcmporaturcs and cho 
nodelcss variables, rcspecslvelyt Element interpolation functions 
N^, 1 "1,4 are the same as for the conventional bilinear four node 

element given in Eq* (6.3), and i -5,8 are the nodeless variable 
Interpolation functions given in Eq. (6.6). 

The combination of the conventional and nodeless variable inter- 
polation functions, Eq. (6.7), provides a quadratic temperature 
distributl<in over the element but with only four element nodes. 
InterelbCir' compatibility is preserved since adjacent elements have 
a common nodeless variable on adjoining edges. The magnitude of the 
nonlinear variation on an element edge is measured by the correspond- 
ing nodeless variable. Temperature distributions for the conventional 
bilinear element and the nsdeless variable element are compared in 
Fig. 18. 


6.? Nodeless Variable Finite Element Formulation 
for Thermal-Structural Analysis 

In this section, the thermal finite element formulation for two- 
dimensional linear transient analysis is described. The formulation 
is valid for bpth the conventional element and the nodeless variable 
element. A four node structural element which will be used in 
junction with the thermal element for computation of thermal stresses 
is also presented. 

6.2.1 Linear Thermal Analysis 

In two-dimensional transient heat conduction, the governing 
differential equation for the temperature distribution T(x,y,t) 
may be expressed in the form of 
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Fig. 18. Two-dimensional element interpolation 
functions. 
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where k and k are the thermal conductivities in x and y 
A y 

directions, respectively, Q is the Internal heat generation rate 
per unit volume, p is the density, and c is the specific heat. 

To derive the element equations and element matrices, the method 
of weighte^iii residuals (see section 2.3) is applied to the governing 
differential equation (6.8). VTith the boundary conditions of 
specified temperatures, surface heating and surface convection as 
shovm in Eqs. (2.3a-c), typical element equations have the form 

[C] {t} + (K^ + K^] (T) - {Q^} + {Qq} + {Q^} + {Qjj} (6.9) 

where [c] is the element capacitance matrix; and are 

element conductance matrices corresponding to conduction and convec- 
tion, respectively. These matrices are expressed in the form of 
integrals over the surface area A of an element with the thickness 
t as follows; 


[ci » t j 

1 pc {N } [N J dx dy 
A 

(6.10a) 

[k„] - t 

C j 

[b ( k] [b„] dx dy 
A '■ 

(6.10b) 

■ [ 

h{N^} [nJ dx dy 

(6.10c) 


where [b,^] denotes the temperature gradient Interpolation matrix, 
is the convection coefficient. The right-hand side of the 


and h 


discretized equation (6.9) contains heat load vectors due to 
specified nodal temperatures, internal heat generation, surface 
heating, and surface convection. These vectors are defined by 


U1 


(q-n){N^}ds (6.11a) 

h 

{Qq}-t| Q{N^}dxdy (6.11b) 

{Q } “ 1 q{N_) dx dy (6.11c) 

q \ i 

(Qh) " j* h {N^} dx dy (6. lid) 


where q is the vector of conduction heat flux across boundary 
that is required to maintain the specified nodal temperatures, q is 
the surface heating rate per unit area, and is the convective 

medium temperature. 

As mentioned earlier, a typical quadrilateral element in 
Cartesian coordinates (x,y) is transformed to the natural coordi- 
nates (5,n) to perform the element matrix integration. In 
computation of the conduction conductance matrix (Eq. (6.10b)), for 
example, the chain rule is first applied to relate the temperature 
gradients in both coordinate systems. 
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Using the coordinate trasnformatlon shown In Eq. (6.2), the above 
relations become 
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where [j] Is the Jacobian matrix defined by 
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(6.14) 


Substituting the element temperature, Eq. (6.1) or (6.7), Into the 
right-hand side of Eq. (6.13) yields 
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(6.15) 


where r is the number of the element unknowns; r = 4 and 8 for 
the conventional bilinear element and the nodeless variable 
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eleraonc, resspectivaly. The temperature gradient interpolation matrix 
in the above equation is given by 


3N, 3N„ 




3n 3n 




(6.16) 


Using dx dy " |j| d^ dn where |j| is the determinant of [j] , 
the conduction conductance matrix terras of the natural coordinates 
is 


C\] - t 


1 1 

I I [k] Ul d? dn 

-1 


(6.17) 


Next, the coefficients in the conduction conductance matrix 
are computed by numerical intecrration; the Lagendre-Gauss method 
is used where the above conduction conductance matrix is written in 
the form, 

NG NG „ 

[K^] " '^i [B^Cq.Hj)] ld(q.nj)l (6.18) 

* 

where W^j^, denote Gauss weight factors, 5^, denote gauss 
integration points and NG is the number of Gauss points in each 
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coordinate direction. Gauss weight factors and Gauss Integration 
point coordinates can be found in Ref. [IS]. 

Other element matrices shown in Eqs. (6.10 -6.11) can be 
formulated in the same manner. For example, the conductance 
matrix and the heat load vector associated with surface convection 
are expressed as 


[Khl 


NG NG 
h 2 2 

i-1 j»l 


5^1 Wj [N^CE^.rij)] lJC5^,nj)l 


(6.19a) 


NG NG 
{Qh> - hT„ E Z 
" i-1 j«l 




(6.19b) 


In performing the numerical integration, the accuracy of the 
matrices depends on the number of Gauss points used. In general, 
the use of n Gauss points provides exact integration when the 
integrand contains polynomials of order up to 2n -1. For the 
conventional bilinear finite element, two Gauss points (NG - 2) 
in each coordinate direction are normally used. Since the nodeless 
temperature interpolation functions contain higher order of 
polynomials than those for the conventional bilinear element, 
more Gauss points should be used. For the linear thermal analysis 
presented herein, three Gauss points in each coordinate direction 
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was found by numerical tests to be appropriate for accurate node- 
less variable element matrices. 

After element matrices are computed, typical element equations 
can be written In the form 

[C] IT) + [k] it) - (Q) (6.20) 

The conventional bilinear element has four nodal temperatures as 
the unknowns, thus the above element equations contain four 
unknowns. The nodeless variable element has four nodal tempera- 
tures and four nodeless variables as the element unknowns, therefore, 
the element equation t contain eight unknowns. In transient 
analysis, these eight equations must be solved simultaneously 
similar to the one-dimensional nodeless variable described in the 
preceding chapter. In steady-state analysis, the four nodelesj 
variable unknowns can be eliminated from the element equations 
using the matrix condensation technique [35], The final number 
of element equations thus reduces to be the same as of the conven- 
tional bilinear element. 

6.2.2 Structural Element 

In this section, the congruent structural element is briefly 
described. The element contains four nodes and permits the same 
discretization with the conventional and nodeless variable thermal 
elements described in the preceding sections. The element stiff- 
ness matrix is the same as used in conventional four node structural 
elements, however, the improved element temperature distributions, 
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Eq. (6.7), are incorporately consistently in the thermal force 
vector computation to yield an integrated thermal structural 
element. 

The structural element at each node has two in-plane displace- 
ment unknowns u and v which may vary with the element local 
coordinates x, y and time t. Element displacement distributions 
are assumed in the form (see Eq. (2.23)), 


» [Ng]{6> 

( 6 . 21 ) 


where Nj, , i *1,4 are the element displacement interpolation 
functions which have the same form as for the conventional finite 
element temperature interpolation functions shown in Eq. (6.3). 

For the quasi-static analysis, typical element equations shown 
in Eq. (2.26) reduce to 

[Kg] m =• {F^} (6.22) 

where [K ] is the element stiffness matrix, and {F } is the 
s T 
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equivalent nodal thermal load vector. These matrices are expressed 
in the form of integrals over the element volume V as 



(6.23a) 


{F^} - J [D] {o}(T(x,y,t) - T^gP dV (6.23b) 


where [B ] is the strain-displacement interpolation matrix obtained 
s 

from the strain-displacement relations, 


'xy 


N. ✓ 


3u 

3x 

3y 

3y 

3u . 3v 
3y 3x 


- [Bg] {6} 


(6.24) 


[d] is the elasticity matrix defined by (plane stress) , 


[D] 


1-v 


1 V 

V 1 

0 0 


0 

1-v 


(6.25) 


where v is Poisson's ratio. The vector {ci} contains the thermal 
expansion coefficients given by (plane stress) 


{a} 


a 
a > 


(6.26) 


0 
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T(x,y,t) is the element temperature computed using the conventional 

or the nodeless variable thermal element and is the reference 

temperature for sero stress. The elasticity matrix [D] and the 

vector of thermal expansion coefficients shown in the above equations 

2 

can be used for plane strain by replacing E/(l-v ) for E, v(l-v) 
for V, and (l+v)a for a. 

Similar to the quadrilateral thermal element, numerical integra- 
tion is required for computing the element matrices. Using the 
Lagendre-Gauss method, the element stiffness matrix and the equivalent 
thermal load vector shown in Eqs. (6.23a-b) are written in the form, 


'■'s' 


NG NG 


E E Wj[Bg(5^,r,j)]" [D] (6.27a) 

i*l j *1 


NG NG T 

{F„} - t E E W. W. [B (I ,nJ]^ 
^ i-1 j«l ^ ^ ® ^ ^ 


[D] {ct} (T(q,nj) lJ(q,nj)l 

(6.27) 


where T(?^,rij) is the temperature at the element Gauss integration 
point 5^ and ri^ . 

Unlike the thermal finite element previously described, the nodal 

• * * 
displacement unknowns of the structural element are the vector 

quantities. Transformation of the element matrices from the local 

coordinates (x,y) to the global coordinates (X,Y,Z) is required, 

t 

In three-dimensions, the element stiffness matrix becomes a 12 by 12 
matrix and similarly with the nodal force vector. Thus the element 
equations contain a total of 12 equations with 12 nodal displacement 
unkonwns in the global coordinates. After the global element matrices 
are assembled and the nodal displacements are computed, element nodal 



dlsplaceaencs In the local coordinates can be obtained. Then the 
element stresses can be computed from 
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/• 


'X 




- {6> - {«} (T(x,y,t) - 


( 6 , 28 ) 


6,3 Applications 


To illustrate the performance of the two-dimensional nodeless 
finite element presented in section 6,1,2, two examples are analyzed: 

(1) a rectangular plate with surface convection, and (2) a simplified 
wing section with aerodynamic heating. In each example, benefits 

of the nodeless variable finite element are demonstrated by comparison 
with results from conventional finite element and analytical solutions. 

6.3.1 Steady-State Heat Conduction in a Plate with Surface 
Convection 

A rectangular plate (Fig. 19(a)) has a specified temperature 
Tq along the boundaries. The plate is cooled by surface convection 
to a zero medium temperature, T„ - 0. Using symmetry, a quarter of 
the plate is first modeled by: (1) one conventional element, and 

(2) one nodeless variable element. 

Fig. 19(b) shows ‘the comparative temperature distributions at 
y « b/2 for an analytical solution [27] , the conventional element 
and the nodeless variable element solutions. For these models, the 
conventional element gives a relatively high error compared to the 
nodeless variable element. The largest error for both finite element 
models occurs at the center of the plate (16% and 3% for the 
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(a) PLATE COOLED BY SURFACE CONVECTION 


Fig. 19. Conventional and nodeless variable finite 
element solutions for a plate with surface 
convection. 
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(b) COMPARATIVE TEMPERATURE DISTRIBUTIONS ALONG y » -g* 


Fig, 19. Conventional and nodeleas variable finite element 
solutions for a plate ^d.th surface convection. 
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conventional element and the nodeless variable element, respectively). 
At the center of the plate, both elements shov a discontinuity of 
conduction heat flux Indicating a need for mesh refinement. Next, 
the plate Is modeled by using four finite elements shown by the dotted 
lines In Fig. 19(a); comparative temperature distributions are shown 
In Fig. 19(c). Four conventional elements provide a fair estimate 
of the temperature variation, but four nodeless variable elements 
yield excellent predictions for both nodal and element temperatures. 
Comparisons of temperatures at other sections of the plate (not shown) 
demonstrate that four nodeleas variable elements provide excellent 
agreement with the analytical solution for the entire plate. 

6.3.2 Simplified Wing Section with Aerodynamic Heating 

To demonstrate the usefulness of the two-dimensional nodeless 
variable elements In aerospace thermal-structural analysis, a simpli- 
fied wing section is analyzed (Fig. 20(a)). Top and bottom skins 
of the wing section are connected by three corrugated spars and are 
subjected to symmetrical, nonuniform time-dependent aerodynamic 
heating. 

Three finite element models are employed to computed temperatures. 
For a unit depth in the spanxd.se direction, the first model consists 
of seven conventional elements; two elements each for the top and 
bottom skins and one element for each spar. The second model is 
identical to the first model except nodeless variable elements are 
used. The third model use's a refined mesh (not shoxm) xdth 35 conven- 


tional elements. 


aRieiWAL PACI&* 

OF POOR QUALITY 


ANALYTICAL 

Ac ▲ CONVENTIONAL F.E., 4 ELEMENTS 

O— -O NODELESS VARIABLE F.E., 4 ELEMENTS 



X. 

■q 


(c) COMPARATIVE TEMPERATURE DISTRIBUTIONS ALONG y » . 


Fig. 19. 


Conventional and nodeless variable finite element 
solutions for a plate with surface convection. 
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(a) SIMPLIFIED WING SECTION WITH AERODYNAMIC HEATING 



Fig. 19. Conventional and nodeless variable finite element 
solutions for a simplified wing section with 
aerodynamic heating. 
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Comparative skin temperature distributions at t « 150 s. are 
shown in Fig. 20(b); the number of elements cltJd is for the skin 
only, The nodeless variable finite element model predicts a realistic 
temperature distribution and gives very good agreement with the 
result from the refined conventional finite element model. The crude 
conventional finite element model underestimates the average skin 
temperature and is unable to provide details of the nonuniform 
temperature distribution. 

in computation of the skin thermal stress, classical beam 
theory [16] is employed for comparison with two finite element stress 
analyses. Detailed temperature distributions from the refined 
conventional finite element thermal model are used to compute the 
stress a from beam theory. Temperature distributions from the 

A 

crude conventional thermal finite element model and the nodeless 
variable thermal finite element model are transferred to a structural 
finite element model with the same discretization for the stress 
computations. Comparative stress distributions at t » 150 s. are 
presented in Fig. 20(c). The advantage of using the improved 
temperature distributions from the nodele.'s variable finite element 
model in computing stresses is clearly demonstrated. These stress 
distributions are in excellent agreement with the result from beam 
theory with both the critical stress and its location accurately 
predicted. Using the temperature distribution from the crude 
conventional finite element model yields significant errors in the 
stress distribution and is unaccetable for this problem. 

These two examples clearly demonstrate the benefits of the 
two-dimensional nodeless variable finite element that can be obtained 
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(b) COMPARATIVE TEMPERATURE DISTRIBUTIONS 
ALONG CHORDWISE DIRECTION, t = ISO 3. 


Fig. 20. Conventional and nodeless variable finite element 
solutions for a simplified wing section with 
aerodynamic heating. 
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(c) COMPARATIVE THERMAL STRESS DISTRIBUTIONS 
ALONG CHORDWISE OIRECTION, t » 150 s. 


Fig. 20. Conventional and nodeless variable finite element 
solutions for a simplified wing section with 
aerodynamic heating. 
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for thermal-structural analysis. Additional applications, a summary 
of the nodeless variable approach and the thermal-structural finite 
element formulation presented in this chapter appear in Ref, [36], 



Chapter 7 


CONCLUDING REMARKS 

An integrated approach for improved thermal-structural finite 
element analysis is presented. The approach was motivated by aero- 
space applications to improve thermal-structural finite element 
analysis capabilities. An important goal is to eliminate the 
incompatibility between thermal-structural analyses where a more 
detailed finite element model is required for the thermal analysis 
than for the structural analysis. The integrated approach is 
characterized by: (1) thermal and structural finite elements 

formulated with common geometric discretization for full compatibility 
during the coupling of the analyses, (2) accurate nodal and element 
temperatures provided by improved thermal finite elements, and (3) 
accurate thermal loads for the structural finite element analysis to 
further Improve accuracy of the structural response. 

Basic concepts and procedures of the integrated thermal-structural 
finite element analysis are described. New thermal finite elements 
for improved thermal analysis accuracy are developed. Thermal finite 
elements which yield exact nodal and element temperatures for one- 
dimensional linear steady-state heat transfer problems are presented. 
These thermal finite elements are formulated based upon using closed- 
form solutions of the governing differential equations. For general- 
heat transfer problems where closed-form solutions are not available, 
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improved thermal finite elements are developed by employing the 
nodeless variable formulation. The nodeless variable finite element 
uses extra unknowns as element variables to permit higher order 
element temperature interpolation functions. Detailed element 
temperature distributions are obtained without using additional 
element nodes while a common discretization with lower order congruent 
structural finite elements are maintained. 

Nodeless variable finite elements are formulated for the 
following heat transfer cases: (1) one-dimensional linear transient 

analysis, (2) one-dimensional nonlinear transient analysis with 
radiation, and (3) two-dimensional linear transient analysis. 

General formulations of the nodeless variable finite elements for 
each heat transfer case are described in detail. For comparison, 
conventional finite elements customarily used in typical finite 
element programs are also presented. Results of temperatures obtained 
from the thermal analysis are transferred directly to the structural 
analysis to compute displacements and stresses. 

To demonstrate the capabilities and efficiency of the Integrated 
finite element approach, several examples in academic and more 
realistic problems are employed. The accuracy of the approach is 
evaluated by comparisons with analytical solutions and conventional 
thermal-structural analyses. Results indicate that the integrated 
finite element approach provides a significant improvement in the 
accuracy and efficiency of thermal-structural analysis and offers 
potential for applications to other coupled problems. 
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APPENDIX A 

EXACT FINITE ELEMENT INTERPOUTION FUNCTIONS 


Exact element interpolation functions in the form of equation 

(3.14) for the thermal finite element Cases 1-8 (Figure 4 and Table 1, 

pp. 39 and 40) are presented. Nodeless parameters are shovm in 

« 

Table 2 (p. 4l) . The lower case letters in parentheses denote heat 
load cases defined in Table 1, General solutions to the differential 
equations for Cases 6 and 7 appear in reference [37]. 


Rod (Case 1) 


N 


1 - 



N. 


X 

I 


N 


0 



M sinh m(L~x) „ _ slnh mx 

" slnh mL ^2 sinh mx 


N 


0 



N 


2 


where m ■ /hp/kA. 


Slab (Case 2) 


N 


1 




N 


0 



(a,c,d) 


(b) 


(a,c) 


* ■I'.'Sm 


16G 


Hollow Cylinder (Caae 3) 


% - 7 »2 - 7 

•*0 “ ^ ^ " 


where w 


ln(i) 

E 


Hollow Sphere (Case 4) 


. a(b~r) 
1 r(b-a) 


N. 


b(r~a) 

r(b-a) 


No ■ ~(r-a) <b-r) (r+a+b) 


Cylindrical Shell (Case 5) 


Ni-1-^ 


^2 " L 


N 


0 L 


fa - f) 


N, - 


slnh m(L-s) 
slnh mL 


N 


_ slnh ms 
2 “ slnh mL 


Nq - 1 - Ni - N2 


(a,c) 

(a,c) 

(h,c) 

(e> c ) d) 
(b) 


where m ■ /h/kt , 
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Conical Shell (Case 6) 


■*1 - 5 


N« " ~ ln(~) 
2 w a 


(a,o,d) 




iQ(ma) KQ(mb) - IgCmb) KqChis) 
^1 “ IgCma) KQ(mb) - iQ(mb) RgCma) 


<b) 


iQ(ma) KQ(ms) - IgCms) RgCma) 
2 " IgCma) KQ(mb) - IqW KQ(ma) 


No ■ 


i - N2 


where w ■ ln(—) , ra ■ /h/lct j Iq and Kg are modified Bessel 
funcClons of the first and second kind of order zero, respectively. 


Spherical Shell (Case 7) 


Nj - 1 - Nj 


In l+sin(s/a) 
l-8in(s/a) 
l+3in(L/a) 
l-sin(L/a) 


i 

(a,c,d) 


Ng ■ In [cos (s/a)] - N2 In [cos(L/a)] 
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Flow Passage (Case 8) 


N 


1 



^2ox 

e 

2aL 

e 


1 

1 


2ax 

e 

^2aL 
- e 


(a,d) 


N 


0 



N 


2 


M . slnh 6(L-x) 
”l ® slnh 8L 


M - slnh. 6x 

”2 slnh 8L 


(b) 


"o ■ I - “l - “2 


where ct “ mc/2kA, 8 ■ and n ■ /hp/kA . 
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APPENDIX B 

FINITE ELEMENT MATRICES FOR ONE-DIMENSIONAL 
LINEAR STEADY-STATE PROBLEMS 

Exact finite element 'oatrlces for the thermal and structural 

» 

finite elements descrlbfiid In Chapter 3 are presented. Theinnal 
conductance matrices and heat, load vectors are given In the form of 
equation (3.27) for Cases 1-6 and equation (3.17) for Case 8 
(Figure 4 and Table 1, pp. 39 and 40 ) . These finite element 
matrices are derived using the exact element interpolation functions 
shown in Appendix A. Similarly, structural stiffness matrices and 
equivalent nodal forces due to thermal loads are derived using the 
element displacement interpolation functions shown in Tables 3 and 
4 (pp. 61 and 67). The lower case letters in parentheses denote 
heat load cases defined in Table ,j, , 
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THEFHAL FINITE ELEMENT MATRICES 

R od (Case 1) 

Conductance Matrices'. 

^11 “ ^22 “ <a,c,d) 

^11 " ^22 " mL)/(m sinh mL) ' (b) 

K ^2 ■ " bp/(m sinh mL) 

Heat Load Vectors: 

- Q 2 - QAL/2 (c) 

“ Q 2 “ qpL/2 (d) 

where m » /hp/kA . 

Slab (Case 2) 

Conductance Matrices: 

^11 " ^22 °* (siC) 


- k/L 
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Heat Load Vector 

- Q 2 - QL/2 (c) 

Hollow Cylinder (Case 3) 

Conductance Matrices 

^2 " ” 

Heat Load Vector 

- Q (- a^/2 + (b^- a^) / 4w) (c) 

Q 2 » Q (b^/2 (b^- a^) / 4w) 

where w = ln(b/a) 

Hollow Sphere (Case 4) 

Conductance Matrices 

^11 “ ^22 ~ kab/(b-a) (a,c) 

K ^2 “ ~ kab/(b-a) 

Heat Load Vector 

= Qa(2a^ + b^ - 3a^b) /(6(&-a)) (c) 
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($2 ■ Qb(a^ + 2b^ - 3ab^) /(6(b-a)) 

Cylindrical Shell (Case 5) 

Conductance Matrices 


hi ■ hi - 

“ ^22 ■ (b cosh mL) / 
Kj ^2 " " b/(mt sinh mL) 
Heat Load Vectors 


(mt sinh mL) 


Ql “ Q 2 “ QL/2 
Qj " Q 2 * qL/2t 

" Q 2 “ hT^ (cosh mL - 1) / (mt sinh mL) 


where m = /h/kt . 


(a,c,d) 


(b) 


(c) 

(d) 
(b) 


Conical Shell (Case 6) 

Conductance Matrices 

Kii = K22 = k/w 



(a,c,d) 


- k/w 
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Heat Load Vectors 

- Q(-a^/2 + (b^-a^)/4w) 

Q 2 - Q(b^/2 - (b^«-a^)/4w) 

■ q(-a^/2t + (b^-a^)/4wt) 
Q 2 - a(b^/2t - (b^-a^)/4wt) 
where w ■ ln(b/a) 


Flow Passage (Case 8) 

Conductance Matrices 

Kg - K- «< kaa(e^“^+ l)/(e^“^ -1) 

11 ‘■22 

Kc - K- « - K_ 

12 ^^21 ‘^11 

=■ kA((BH/ZG) - (a/2) - (3^E(E+F)/2aG^)) 

K» =» kA((3H/2G) + (a/2) - (3^E(E-F)/2aG^)) 
22 

Kc = Kc = kA(-(3^EH/2aG^) - (3F/2G)) 

12 21 

K-j “ =» - me/ 2 

'^ll ^12 


K., 


21 


K = mc/2 
^22 


( 0 ) 


<d) 


(a,d) 


(b) 


(a.d) 
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fi 


K„ - - K - - mc/2 (b) 

^11 ^^22 

K„ - - - gE/oG 

^12 21 

K. - hp((pH/G) + (a) - (e^E<F+E)/aG^)) /2(g^ -a^) (b) 

K. - hp<(0H/G) - (a) - (E^E(F-E)/aG^)) /2(E^ -a^) 

^22 

K. - Kt, - hp((sWaG^) - (6F/G)) /2(0^ -a^) 

"21 ■ 

Heat Load Vectors 

Qj^ - qp(l -e^“^+ 2«L) /2a(l (d) 

f^L / i ^ ^cxL rs • 2otLv i M 2ccLv 

® qp(-l -he -2aL e ) /2c*^l -e ) 

- hpT„(e(H+E-F) -aG) /G(p^ -a^) (b) 

Q 2 =• hpT^CE(H-E-F) +aG) /G(g^ -a^) 

/~2 T 

where a ** ic/2kA, $ * f a + m , m =• vhp/kA , E “ sinh aL , 


F =» cosh aL , G » sinh 3L , H “ cosh BL . 
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STRUCTURAL FINITE ELEMENT MATRICES 

Truss Element (Case 1) 

Stiffness Matrices 

*^11 “ ^22 " (a,b,c,d) 

Ki 2 - K22 ■ - 

Force Vectors 

Fj_ - - F 2 » - aEA(TQ/6 + (Tj_+T 2 )/ 2 ) (a,c,d) 

Fi - - F 2 - - aEA(Cj_Tg + C2(Tj^+T2)) (b) 

where Cj_ » 1 - (2(cosh mL - l)/mL sinh mL) . 

C 2 * (cosh mL - l)/mL sinh mL 

Axisymmetric Element (Case 3, Plane Stress) 

Stiffness Matrices 

= E((b^+a^) - v(b^- a^))/(l-v^)(b^'-a^) (a,c) 

K22 - E((b^-a^) + v(b^-a^))/(l-v^)(b^-a^) 

Kj _2 = E(>2ab)/(1 - v^)(b^ -a^) 



Force 


where 
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Vectors 

- - aEaP/2w(l -v)(b^-a^) (a,c) 

- bEaP/2w(l-v)(b^-a^) 

P - ( -2a^w + b^ - a^)(Tj^ + a^'w T^/ b^) 

+ (2b^w - b^ + a^)(T 2 + w Tq) 

- (b^ - aS Tq / b^ - 2(b^- a^) wT^gf 
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APPENDIX C 

FINITE ELEMENT MATRICES FOR ONE-DIMENSIONAL 
LINEAR TRANSIENT PROBLEMS 

Finite element capacitance matrices for the thermal rod and 
axisymmetric elements described in Chapter 4 are presented. The 
conductance matrix coefficients Kqq, and heat load vector 
components are presentad; the coefficients and Q^, 

i,j * 1,2 appear in Appendix B, Cases 1 and 3. The lower case 
letters in the parentheses denote heat load cases defined in Table 1, 
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Rod Element 

Capaclcance Matrices 
Cqq - pcAL/30 

Cqi « Cq 2 “ pcAL/12 
°11 " ^2 “ 

/ 

“ pcAL/6 

Coo " pcA(((cosh mL - l)/sinh mL)(L/8±nh mL - 3/m) + L) 

2 

^01 " ^02 * P<^A.((1 - cosh mL) (mL - sinh mL)/2a sinh mL) 

2 

^11 " ^22 " PcA((sinh mL cosh mL - mL)/2m sinh mL) 

2 

Ci 2 » pcA((mL cosh mL - sinh mL)/2m sinh mL) 

Conductance Matrices 

Koo - W3L 

>^00 - (h?/m)(mL - 2(cosh mL - l)/sinh mL) 

Heat Load Vectors 

Qq “ QAL/6 
Qq ” qpL/6 


(a^ C| d) 


(b) 


a, c I d) 

(b) 

(c) 

(d) 


Qq ■ hpT„ (L - 2 (cosh mL - l)/m slnh roL) 


Vr> 

(b) 

, where m ■ /hp/tcA 

Axlsymmetrlc Element 

Capacitance Matrices 

Cqq - pc(b^-a^)(4w^(aV a^b^ + b^) + 9w(aW) (a,c) 

+ 6(aW)^)/24b^ 

- - pc(4aV + w(7a^+3b^)(a^-b^) + 4(a^-b^)^)/16wb^ 

Cq 2 = pc(4bV w(7b^’{=33^)(a^-=b^) + 4(a^-b^)^)/l6wb^ 

^11 “ Pc(b^-a^(l+2w+2w^))/4w^ 

^12 " Pc(a^"b‘" + w(a^-b^))/4w^ 

^22 " Pc(b^(l-2w+2w^) - a^)/4w^ 

Conductance Matrices 

Kqq - kw(w(l-(a/b)'^) - (l-(a/b)^)^) (a,c) 

Heat Load Vector 

Qq - (Qb^/4) (w(l-a^/b^) - a-B.^/hV) 


(o) 
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APPENDIX D 

FINITE ELEMENT MATRICES FOR ONE-DIMENSIONAL NONLINEAR 
TRANSIENT ANALYSIS WITH RADIATION HEAT TRANSFER 

Jacobian matrices and residual heat load vectors for the 
conventional finite element and the nodeless variable finite 
element described In Chapter 5 are presented. These element 
matrices which appear In Eq. (5.34) are given In the form of computer 
subroutines. The subroutines are written in FORTRAN IV where the 
definitions of variables used are provided. 


VARIABLES USED IN THE FOLLOWING SUBROUTINES ARE DEFINED 
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ALL RESIDUAL LOAD VECTORS ARE REPRESENTED BY VARIABLES BEGIN 

WITH R f_I E.G. RRAD REPRESENTS RESIDUAL HEAT LOAD VECTOR 

FROM CONDUCTANCE RADIATION MATRIX. 
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